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abstract! 


A  new  compressible  Navier-Stokes  code  in  cylindrical  coordinates  was 
developed  for  investigating  axisymmetric  wakes  of  bluff-based  bodies  in  su¬ 
personic  flows.  In  this  code,  high-order  compact  finite  differences  derived 
for  non-equidistant  grids  are  employed  and  a  new  state-of-the-art  axis  treat¬ 
ment  is  incorporated.  Additionally,  the  fully  three-dimensional  transport 
equations  for  turbulent  kinetic  energy  and  turbulent  dissipation  are  imple¬ 
mented  to  enable  (steady  or  unsteady)  Reynolds  Averaged  Navier  Stokes 
(RANS)  simulations.  Furthermore,  a  new  ’’Flow  Simulation  Methodology” 
(FSM)  was  developed  for  computing  complex  compressible  flows.  The  cen¬ 
terpiece  of  FSM  is  a  strategy  to  provide  the  proper  amount  of  modeling  of 
the  subgrid  scales.  This  is  accomplished  by  a  ’’contribution  function”  which 
locally  and  instantaneously  compares  the  smallest  relevant  scales  to  the  lo¬ 
cal  grid  size.  The  contribution  function  is  designed  such  that  it  provides  no 
modeling  if  the  computation  is  locally  well  resolved  so  that  the  computation 
approaches  a  Direct  Numerical  Simulation  (DNS)  in  the  fine  grid  limit,  or 
provides  modeling  of  all  scales  in  the  coarse  grid  limit  and  thus  approaches 
an  unsteady  RANS  (URANS)  calculation.  In  between  these  resolution  limits, 
the  contribution  function  adjusts  the  necessary  modeling  for  the  unresolved 
scales  while  the  larger  (resolved)  scales  are  computed  as  in  traditional  Large 
Eddy  Simulations  (LES).  Preliminary  results  have  shown  that  the  new  high- 
order  code  has  great  advantages  for  supersonic  base  flow  simulations  and 
that  calculations,  in  particular  together  with  FSM,  will  allow  simulations  of 
supersonic  baseflows  at  much  higher  Reynolds  numbers  than  possible  with 
conventional  LES. 

^Upon  recommendation  from  ARO,  the  scope  of  work  was  shifted  from  passive  and 
active  control  of  supersonic  axisymmetric  base  flows  towards  development  of  a  compact 
Navier-Stokes  code  in  combination  with  a  new  Flow  Simulation  Methodology. 
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1.  INTRODUCTION 

For  axisymmetric  aerodynamic  bodies  in  supersonic  flight,  the  flow  field 
in  the  wake  region  has  considerable  effect  on  the  aerodynamic  drag.  Even 
small  changes  in  the  flow  behavior  of  the  wake  may  affect  the  performance 
of  the  entire  flight  vehicle,  e.g.,  missiles,  rockets,  or  projectiles.  The  effect 
of  the  wake  on  the  aerodynamic  drag  is  mainly  due  to  the  recirculation 
region  that  develops  in  the  base  region  of  the  body  and  thus  to  the  low 
pressure  associated  with  the  recirculating  flow  (base  drag).  Flight  tests  with 
projectiles  (U.S.  Army  549  projectile)  have  shown  that  the  base  drag  may 
account  for  up  to  35%  of  the  total  drag  (see  Rollstin,  1987).  This  suggests 
that  attempts  to  modify  the  near  wake  flow  such  that  the  base  pressure  would 
increase  could  be  highly  rewarding  with  respect  to  drag  reduction  and,  as  a 
consequence,  with  regard  to  increasing  the  performance  characteristics  of 
flight  vehicles  or  projectiles.  However,  in  order  to  modify  the  near- wake  flow 
effectively,  a  detailed  understanding  of  the  underlying  physical  mechanisms 
that  are  responsible  for  the  flow  behavior  is  required. 

For  this  reason,  in  the  past,  numerous  research  efforts,  experimental, 
theoretical,  and  computational,  have  focused  on  understanding  the  physics 
of  the  base  flow  or  the  near-wake  region.  However,  due  to  the  difficulty  of 
the  problem  and  limitations  of  experimental  and  numerical  techniques,  these 
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studies  focused  almost  exclusively  on  the  mean  flow  behavior  and  not  on 
the  potential  effects  of  the  unsteady  structures  of  the  turbulent  flow.  As 
discussed  below,  if  energetic  (coherent)  structures  are  indeed  present  in  such 
flows,  and  more  recent  evidence  indicates  that  they  are  [see  Smith  Sc  Dutton 
(1996),  Bourdon  Sc  Dutton  (1998),  Harris  Sc  Fasel  (1996),  and  Tourbier  Sc 
Fasel  (1994)],  the  dynamical  behavior  of  such  structures  may  have  a  profound 
effect  on  the  resulting  mean  flows.  Also,  as  a  consequence  of  this,  profound 
modifications  of  the  base  flow  behavior  and  thus  the  base  drag  must  involve 
a  modification  of  the  dynamics  of  these  large  structures. 

In  recent  years,  techniques  for  reducing  the  base  drag  have  been  suggested 
and  investigated,  such  as  boattailing,  base  bleed,  base  burning,  etc.  [see,  for 
example,  Sahu  et  al.  (1985),  Nietubicz  Sc  Gibeling  (1993),  Sahu  Sc  Heavey 
(1995),  Sahu  Sc  Nietubicz  (1994),  Mathur  Sc  Dutton  (1996a, 6)  and  Bourdon 
Sc  Dutton  (2000a)].  Such  techniques  can  indeed  modify  the  mean  base  flow 
and,  as  a  consequence,  the  base  drag.  However,  the  underlying  mechanisms 
responsible  for  the  changes  are  not  understood.  In  order  to  optimize  such 
techniques  for  practical  applications,  it  is  essential  that  the  fundamental 
mechanisms  that  are  responsible  for  these  changes  are  brought  to  light.  As 
discussed  below,  the  key  to  uncovering  these  essential  mechanisms  is  in  first 
understanding  the  role  of  the  (coherent)  turbulent  flow  structures  and,  in 
particular,  how  the  dynamics  of  these  structures  is  influenced  by  proposed 
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techniques  for  drag  reduction. 

1.1  Experimental  Research 

Historically,  the  fundamentals  of  the  near-wake  behavior  were  first  inves¬ 
tigated  using  experiments  almost  exclusively.  However,  wind  tunnel  experi¬ 
ments  suffered  from  the  difficulty  of  properly  supporting  the  (axisymmetric) 
base  model  so  as  not  to  cause  undue  effects  on  the  flow  field  behind  the 
base  Chapman  (1951).  In  experimental  investigations,  different  methods  of 
model  support  have  been  proposed;  for  example,  rear  sting  support  Donald¬ 
son  (1955),  side  mounted  struts  Chapman  (1951),  and  up  stream  wire  mounts 
Dayman  (1963).  However,  all  of  these  techniques  have  shown  non  negligible 
effects  on  the  near  wake  behavior  and,  in  particular,  on  the  base  pressure. 
The  forward  sting  support  used  by  Herrin  &:  Dutton  (1991)  and  also  in  sub¬ 
sequent  and  current  experimental  investigations  appears  to  be  superior  to 
other  methods  with  regard  to  minimizing  flow  interference  [also  see  Dutton 
&  Addy  (1993),  Herrin  &  Dutton  (1995),  Herrin  &;  Dutton  (1994),  Mathur 
&  Dutton  (1996a),  Mathur  &  Dutton  (19966)  and  Bourdon  &  Dutton  (1998, 
2000o)]. 

In  addition  to  the  interference  caused  by  model  support,  wind  tunnel 
interference  can  have  a  considerable  influence  on  the  base  flow  behavior. 


Therefore,  extra  care  has  to  be  taken  in  order  to  understand  and  minimize  the 
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effects  of  vibrations,  upstream  sound,  noise  radiated  from  turbulent  boundary 
layers  of  the  wind  tunnel  walls,  blockage,  etc.  This  is  particularly  important 
for  investigating  the  role  of  large  coherent  structures  and  for  reliably  testing 
passive  control  techniques  (see  section  2.2).  For  the  wake  flow  with  a  central 
jet  (power  on  mode),  wind  tunnel  interference  would  be  somewhat  less  of 
a  problem  and  therefore  appears  to  be  experimentally  easier.  However,  for 
this  case,  the  shear  layers  become  much  thinner,  making  it  much  harder  to 
resolve  the  resulting  large  gradients  when  taking  measurements. 

Furthermore,  the  quality  of  experimentally  obtained  flow  data  in  the  wake 
is  often  compromised  whenever  intrusive  probes  (such  as  Pitot  pressure  tubes, 
hot  wire  probes,  etc.)  are  used.  This  is  particularly  true  for  measurements 
within  the  recirculating  region  immediately  downstream  of  the  blunt  base 
of  the  body.  The  strong  sensitivity  to  intrusive  probes  is  due  to  the  ab¬ 
solute/global  instability  that  may  be  present  in  supersonic  wake  flows  (see 
section  2.2).  Absolute  stability  behavior  is  very  sensitive  to  small  changes  of 
the  mean  flow  velocity  proflles,  in  particular  when  solid  surfaces  (probes)  are 
inserted,  which  can  change  the  upstream  pressure  feedback  mechanism. 

In  free-flight  experiments,  many  of  the  above  mentioned  difficulties  do  not 
arise.  However,  in  addition  to  the  much  greater  cost,  free  flight  experiments 
also  suffer  from  major  disadvantages.  For  example,  it  is  not  possible  to  map 
out  details  of  the  flow  field  as  is  possible  with  wind  tunnel  experiments. 


Also,  the  reliability  of  flight  test  data  is  often  in  question  because  of  the 
considerable  difficulties  in  controlling  conditions  and  parameters  in  the  free 
flight  tests. 

In  recent  years,  many  experimental  efforts  have  focused  on  modifying  the 
base  flow  in  order  to  reduce  the  base  drag.  So  far,  only  so-called  passive  con¬ 
trol  methods  have  been  considered,  such  as  boattailing,  base  bleed,  and  base 
burning  [see,  for  example,  Cortwright  &  Schroeder  (1951),  Reid  Sz  Hastings 
(1959),  Bowman  &  Clayden  (1968),  Clayden  &  Bowman  (1968),  Valentine 
&:  Przirembel  (1970),  Hubbartt  et  al.  (1981),  Ding  et  al.  (1992),  Mathur  &: 
Dutton  (1996a, i»)].  These  experiments  have  shown  that  there  are  potentially 
considerable  rewards  for  passive  flow  control.  However,  it  is  still  not  under¬ 
stood  why  certain  measures  are  more  effective  than  others  (or  why  others  do 
not  work  at  all)  and  what  the  optimal  parameters  should  be.  The  reason  for 
this  is  that  the  fundamental  physical  mechanisms  in  a  supersonic  wake  flow 
are  not  yet  understood. 

1.2  Numerical  Simulations 

Because  of  the  difficulties  and  uncertainties  associated  with  experiments, 
researchers  have  looked  to  numerical  simulations  as  attractive  alternatives 
and/or  complements  to  experimental  investigations.  Sparked  by  the  rapid 
increase  in  the  computing  power  of  supercomputers  in  the  last  two  decades. 


numerous  attempts  were  made  to  calculate  the  flow  fleld  around  and  be¬ 
hind  axisymmetric  bodies  aligned  with  the  free  stream.  In  numerical  simu¬ 
lations,  problems  associated  with  wind  tunnel  interference,  model  support, 
probe  intrusion,  etc.,  are  not  present.  In  addition  to  possibly  providing 
further  understanding  of  the  relevant  physics  involved,  these  calculations 
were  motivated  by  the  considerable  challenge  that  this  complicated  flow  field 
provides  for  computational  fluid  dynamicists.  The  computational  challenge 
arises  mainly  from  the  combination  of  shock  waves,  thin  free  shear  layers, 
boundary  layers,  and  recirculating  regions;  associated  with  this  are  highly 
disparate  length  scales  and  local  regions  of  very  high  gradients.  Reliable  and 
realistic  computations  can  therefore  only  be  performed  if  the  high  gradients 
can  be  adequately  resolved. 

Furthermore,  the  wake  flow  is  turbulent  (or  at  least  transitional),  requir¬ 
ing  adequate  turbulence  models.  In  practically  all  earlier  and  in  most  of  the 
more  recent  numerical  attempts,  only  the  steady  flow  fleld  has  been  calcu¬ 
lated.  In  these  computations,  ’’Reynolds  Averaged  Navier  Stokes”  (RANS) 
formulations  were  used  in  combination  with  a  turbulence  model,  such  as 
the  algebraic  eddy  viscosity  model  of  Baldwin  and  Lomax,  the  two  equation 
K  —  e  model,  or  the  K  —  u  model.  In  the  literature,  numerous  attempts 
have  been  reported  on  calculating  turbulent  wake  flows  using  either  the  com¬ 
plete  Navier  Stokes  equations  or  the  thin  layer  (parabolized)  Navier  Stokes 
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equations,  or  combinations  thereof  [for  example,  the  parabolized  equations 
in  part  of  the  flow  fleld  (typically  up  to  the  trailing  edge  of  the  body)  and  the 
full  equations  in  another  part  (typically  in  the  wake  region)].  Also,  variable 
(adaptive)  grids  were  used  in  the  computations  in  order  to  resolve  the  local 
large  gradients  of  the  flow  fleld. 

Relatively  few  direct  comparisons  of  numerical  simulations  with  labora¬ 
tory  experiments  were  attempted  in  the  past.  The  more  successful  numerical 
calculations  were  able  to  capture  the  qualitative  features  of  the  experimen¬ 
tally  obtained  flow  fleld.  However,  relevant  quantities,  such  as  turbulent 
shear  stresses  and,  in  particular,  such  fundamental  aspects  as  the  all  impor¬ 
tant  base  pressure  or  the  size  of  the  separation  bubble  (recirculation  region), 
were  often  not  predicted  well  [see  Petrie  &;  Walker  (1985)].  Also,  and  partic¬ 
ularly  disturbing,  using  different  turbulence  models  but  otherwise  identical 
codes  yielded  inconsistent  (and  often  contradictory)  results.  This  is  an  indi¬ 
cation  that  something  fundamental  might  be  amiss. 

Childs  &  Caruso  (1987)  surveyed  existing  computational  Navier  Stokes 
methods  for  calculating  turbulent  compressible  wake  flows.  They  suggested 
that  inadequate  numerical  resolution  of  the  large  gradients  and,  in  particu¬ 
lar,  the  turbulence  modeling  may  be  the  major  culprits  for  not  being  able 
to  reliably  and  accurately  calculate  the  flow  field.  It  appears  that  because 
of  the  strong  effect  of  different  turbulence  models  on  the  results,  even  the 
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qualitative  agreement  between  various  computations  and  experiments  may 
often  have  been  entirely  fortuitous.  Similar  observations  were  made  by  Sahu 
&:  Nietubicz  (1994),  who  investigated  the  influence  of  different  turbulence 
models  for  RANS  calculations  of  supersonic  axisymmetric  base  flows.  They 
showed  that,  while  some  flow  quantities  can  be  modeled  adequately  by  ad¬ 
justing  certain  parameters  in  the  turbulence  model,  other  quantities  do  not 
agree  well  with  experimental  results. 

Although  considerable  progress  has  been  made  in  past  years  in  using 
RANS  calculations  for  supersonic  axisymmetric  wake  flows  behind  missile 
type  bodies,  there  is  still  no  true  reconciliation  between  experimental  and 
computational  results  [see  Delery  &  Wagner  (1990)].  In  addition  to  possible 
shortcomings  of  the  numerical  codes,  the  discrepancies  were  due  in  part  to  a 
lack  of  reliable  experimental  data,  in  particular  due  to  a  lack  of  experimental 
data  in  sufficient  detail  to  allow  scrutinizing  comparisons  between  experi¬ 
ments  and  calculations.  Now,  due  to  the  extensive  experimental  effort  at  the 
University  of  Illinois  by  Dutton  and  co-workers,  sufficient  and  reliable  exper¬ 
imental  data  are  being  made  available  for  close  comparison  with  numerical 
computations. 

In  spite  of  the  discussed  shortcomings  of  the  Reynolds-averaged  calcula¬ 
tions,  recent  applications  have  attacked  increasingly  more  difficult  situations. 
As  modification  of  wake  flows  is  now  being  considered  as  a  means  of  drag 
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reduction,  RANS  calculations  have  been  performed  to  investigate,  for  exam¬ 
ple,  the  effects  of  base  bleed  [see,  for  example,  Sahu  et  al.  (1985),  Sahu  & 
Heavey  (1995)  and  Danberg  &  Nietubicz  (1992)]  and  base  burning  Nietubicz 
&  Gibeling  (1993).  The  simulations  show  a  correct  trend  for  the  drag  reduc¬ 
tion  when  compared  with  experimental  results.  However,  in  some  cases,  even 
global  flow  field  characteristics  such  as  the  centerline  velocity  are  not  in  good 
agreement  with  experiments.  This  is  an  indication  that  certain  mechanisms 
are  difficult  to  model  in  RANS  calculations.  Prom  the  experimental  effort 
of  Dutton  and  co-workers,  it  is  becoming  increasingly  more  obvious  that  at¬ 
tempts  to  only  calculate  the  steady  flow  might  miss  the  essential  physical 
mechanisms  that  are  responsible  for  the  base  flow  behavior.  Dutton  and 
co-workers  have  identified  a  significant  time-dependent  behavior  and  found 
large  coherent  structures  that  may  play  crucial  roles  in  the  base  drag.  These 
findings  may  be  particularly  important  when  modifications  of  the  flow  are 
attempted  (for  example,  using  boattailing,  base  bleed,  etc.).  Therefore,  in 
order  for  numerical  simulations  to  be  useful  for  uncovering  mechanisms  rel¬ 
evant  to  the  base  drag  in  supersonic  flows,  they  have  to  be  able  to  capture 
the  unsteady  behavior  resulting  from  the  dynamics  of  the  large  (coherent) 
turbulent  flow  structures  (see  discussion  section  2.2). 

Another  approach  of  turbulence  modelling  are  Large  Eddy  Simulations 
(LES).  For  LES,  only  the  large  scale  motion  is  computed  directly.  The  sub- 
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grid  scale  structures  which  are  assumed  to  be  isotropic  are  modelled.  This 
method  was  initially  proposed  by  Smagorinsky  (1963)  who  based  his  model 
on  the  Boussinesq  assumption  (Boussinesq  (1877))  which  states  that  the 
turbulent  stresses  can  be  related  linearly  to  the  mean  velocity  gradients  by 
an  eddy  viscosity. 

For  the  flow  under  consideration,  LES  would  be  the  method  of  choice  be¬ 
cause  of  the  ability  to  capture  three-dimensional  and  unsteady  effects.  The 
shortcomings  of  the  standard  Smagorinsky  LES  model  which  include  the 
failure  to  account  for  normal  stresses  (which  are  present  in  the  recirculation 
region  for  example),  an  ad-hoc  constant  that  is  only  valid  for  a  certain  flow 
region  and  "damping  functions  that  are  needed  at  the  walls  which  are  em¬ 
pirical  in  nature  and  do  not  apply  to  general  wall  bounded  turbulent  flows” 
(Speziale  (1997))  make  this  model  ill-suited  for  this  research.  Further  short¬ 
comings  of  this  widely  used  model  can  also  be  found  in  Jimenez  &  Moser 
(1998). 

In  the  last  few  years,  a  new  class  of  hybrid  RANS/LES  turbulence  mod¬ 
els  has  emerged  that  tries  to  unify  the  best  of  both  worlds.  The  underlying 
principle  is  to  conduct  a  RANS  type  simulation  close  to  a  wall  and  perform 
an  LES  away  from  the  wall  to  capture  the  unsteadiness.  One  way  of  imple¬ 
menting  this  method  has  been  proposed  by  Spalart  et  al.  (1997)  and  is  known 
as  Detached  Eddy  Simulation  (DES).  In  principle,  the  method  compares  the 
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distance  to  the  closest  wall  with  the  product  of  the  local  grid  spacing  and  an 
empirical  constant.  For  small  distances,  the  eddy  viscosity  is  computed  ac¬ 
cording  to  a  one-equation  RANS  model,  otherwise  a  Smagorinsky  type  LES 
is  computed. 

This  model  has  been  applied  to  the  experimental  case  of  Dutton  and 
co-workers  by  Forsythe  &  Hoffmann  (2000)  who  mainly  explored  different 
choices  for  the  empirical  constant.  The  mean  flow  quantities  in  the  recircu¬ 
lation  region  were  found  to  deviate  significantly  from  the  experiments,  even 
though  they  reported  very  good  agreement  of  the  mean  base  pressure  distri¬ 
bution.  Later  work  by  Forsythe  et  al.  (2002)  incorporated  the  compressible 
shear  layer  corrections  into  the  RANS-part  of  the  simulation  and  it  was  shown 
that  this  modification  had  a  negligible  impact  on  the  solution. 

Another  methodology  is  based  on  using  the  same  turbulence  model  in 
both  limits  and  smoothly  blend  from  an  LES  to  a  RANS.  This  is  achieved  by 
a  contribution  function  which  locally  determines  how  much  of  the  turbulence 
is  resolved  by  the  computational  grid.  This  new  Flow  Simulation  Methodol¬ 
ogy  (FSM)  was  proposed  by  Speziale  and  Fasel  (see  Speziale  (1998),  Sandberg 
(1999),  Zhang  et  al.  (2000)  and  recently  extended  to  compressible  flows  by 
Terzi  &  Fasel  (2002)).  The  turbulence  model  is  a  state-of-the-art  explicit  al¬ 
gebraic  Reynolds  stress  model  (ASM)  developed  by  Gatski  &  Speziale  (1993), 
which  is  derived  from  the  full  Reynolds  transport  equations.  It  features  an 


anisotropic  eddy  viscosity  with  strain  dependent  coefficients  and  accounts  for 
normal  stresses  and  rotational  effects. 

The  contribution  function  is  determined  by  comparing  the  local  grid  size 
with  a  local  estimate  of  an  appropriate  turbulent  length  scale,  usually  chosen 
to  be  the  Kolmogorov  length  scale  (Lk)-  For  very  coarse  grids  (compared 
to  Lfc)i  the  Reynolds  stress  from  the  turbulence  model  is  scaled  with  unity 
and  a  RANS  simulation  is  conducted  (Note  that  the  ASM  is  capable  of 
capturing  unsteadiness,  therefore  an  unsteady  RANS  (URANS)  simulation 
is  performed).  In  the  limit  of  fine  resolution,  the  Reynolds  stress  is  multiplied 
with  zero,  recovering  a  DNS.  Between  these  two  limits,  a  ’’non-traditional” 
LES  based  on  the  nonlinear  ASM  is  performed.  This  new  methodology  seems 
to  be  very  promising  for  the  fiow  under  investigation  as  it  will  enable  to 
capture  the  relevant  large  structures  and  should  therefore  be  able  to  predict 
the  mean  fiow  quantities  in  the  near- wake  region. 

Towards  this  end,  we  initiated  several  years  ago  (with  funding  from  ARO) 
a  computational  program  that  exactly  addresses  these  issues  [see  Tourbier 
&  Fasel  (1994),  Tourbier  (1996),  and  Harris  &  Fasel  (1996)].  We  believe 
that  high-quality  Direct  Numerical  Simulations  (DNS),  Large- Eddy  Simu¬ 
lations  (LES),  unsteady  RANS  simulations  (URANS)  and  especially  FSM 
are  required  to  uncover  the  essential  physics  and,  equally  important,  due  to 
the  complexity  and  difficulties  in  experiments  and  simulations  as  discussed 
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above,  that  these  simulations  are  carried  out  simultaneously  and  in  close 
collaboration  with  experiments. 
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2.  TECHNICAL  BACKGROUND 

2.1  Mean  Flow  Field  Behind  a  Cylindrical  Body  with  a  Blunt 
Base 

The  mean  flow  field  for  the  near  wake  region  downstream  of  an  axisym- 
metric  body  with  a  blunt  base  is  shown  schematically  in  figure  2.1.  This  flow 
topology  is  based  on  the  experiments  by  Dutton  and  co-workers  Herrin  &; 
Dutton  (1994).  This  mean  flow  field,  shown  in  figure  2.1,  disguises  the  fact 
that  the  flow  is  turbulent  and  highly  unsteady.  However,  even  the  mean- 
flow  topology  is  already  very  complex:  Downstream  of  the  expansion  waves 
emanating  from  the  sharp  corner,  a  thin,  yet  very  energetic  shear  layer  de¬ 
velops  that  encloses  a  massively  separated  (recirculation)  region  with  a  rear 
stagnation  point.  The  shear  layer  divides  the  flow  field  into  a  high-speed 
supersonic  outer  region  and  a  low-speed  region  near  the  base.  Downstream 
of  the  rear  stagnation  point,  recompression  occurs  and  the  flow  field  from 
there  on  looks  increasingly  like  that  of  a  subsonic  wake  at  several  diameters 
from  the  base.  The  behavior  in  the  separated  (recirculation)  region  has  a 
considerable  influence  on  the  base  drag,  yet  this  region  is  strongly  dependent 
on  the  topology  of  the  shear  layer  and  the  ensuing  recompression  region.  The 
sharp  radial  gradients  in  the  shear  layer  and  the  rapid  spatial  changes  make 
the  flow  field  highly  complex,  and  it  is  therefore  very  difficult  to  investigate 


both  experimentally  and  computationally,  even  if  the  flow  field  were  steady. 
However,  as  mentioned,  the  mean  flow  topology  of  Figure  2.1  is  misleading 
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Figure  2.1  Schematic  of  flow  field  behind  a  cylindrical  body  with  a  blunt 
base  (mean  flow). 

and  disguises  the  fact  that,  in  reality,  the  flow  is  highly  unsteady,  caused  by 
the  very  energetic  shear  layers,  the  flow  recirculation,  and  rapid  downstream 
adjustment  in  the  recompression  region. 

Dutton  and  co-workers  have  started  to  document  the  unsteady  nature 
of  this  flow  and  have  found  that,  even  in  the  dead-air  region,  significant 
dynamical  behavior  of  the  flow  can  be  observed  with  velocities  up  to  100 
m/s  (Herrin  &  Dutton  (1994)).  The  strong  dynamical  behavior  in  the  (mean) 
recirculating  region  was  also  studied  in  detail  using  our  computational  tools 
[Tourbier  &  Fasel  (1994),  Harris  &  Fasel  (1996);  see  also  section  6].  The 
question  is,  of  course,  how  relevant  the  unsteady  nature  of  the  flow  is  to  the 
performance  of  projectiles  and  missiles. 


24 


2.2  Large-Scale  Coherent  Structures 

It  is  well  known  that  for  subsonic  (incompressible)  wakes,  the  dynamics 
of  the  large  (coherent)  structures  play  a  dominant  role  in  the  local  and  global 
behavior  of  the  flow.  This  evidence  was  found  from  both  experimental  inves¬ 
tigations  and  numerical  simulations  (including  ours)  and  was  confirmed  by 
theoretical  studies.  For  incompressible  bluff  body  wakes,  it  has  been  well  es¬ 
tablished  that  the  existence  of  absolute  and  global  instabilities  is  responsible 
for  the  development  of  the  large  structures  (Huerre  &  Monkewitz  (1990)). 
Using  numerical  simulations,  absolute/global  instabilities  were  found  for  a 
two-dimensional  bluff  body  with  a  blunt  base  by  Hannemann  &  Oertel  (1989) 
and  for  an  axisymmetric  body  with  a  blunt  base  by  Schwarz  et  al.  (1994). 
The  absolutely/globally  unstable  modes  for  the  axisymmetric  base  flow  are 
of  a  helical  nature. 

For  supersonic  speeds,  on  the  other  hand,  relatively  little  is  known  about 
the  dynamical  behavior  of  the  large  structures  in  turbulent  flows  or,  in  par¬ 
ticular,  if  absolute/global  instabilities  exist.  This  is  true  for  supersonic  flows 
in  general  and  for  axisymmetric  wakes  in  particular.  The  explanation  for 
this  void  is  that  experiments  are  difficult  to  conduct  for  supersonic  speeds 
(see  discussion  in  introduction).  In  addition,  expensive  facilities  and  intri¬ 
cate  diagnostic  equipment  are  required  for  supersonic  hydrodynamic  stability 
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experiments. 

Due  to  the  lack  of  guidance  from  experimental  investigations  prior  to 
our  computational  research  initiative  on  the  subject,  no  successful  computa¬ 
tional/theoretical  attempts  have  been  made  to  study  the  unsteady,  dynami¬ 
cal  behavior  of  transitional  or  turbulent  supersonic  axisymmetric  base  flows. 
We  feel  that,  as  for  subsonic  wakes,  it  is  exactly  the  dynamical  behavior  of 
the  large  scale  structures  that  strongly  influences  the  flow  fleld  in  supersonic 
axisymmetric  wakes.  We  conjecture  that  the  discrepancies  between  experi¬ 
ments  and  steady  Reynolds  averaged  Navier  Stokes  calculations  may  be  due 
to  the  fact  that  this  unsteady,  dynamic  behavior  of  the  large  scale  structures 
is  not  properly  accounted  for  in  current  turbulence  models.  The  fact  that 
the  results  of  Reynolds  averaged  calculations  depend  strongly  on  the  turbu¬ 
lence  models  used  [see  Petrie  &:  Walker  (1985)  and  Sahu  (1992)]  supports 
this  conjecture. 

As  in  the  subsonic  case,  supersonic  wake  flows  behind  a  blunt  base  are 
dominated  by  a  separated  region,  and  associated  with  it  is  (in  the  mean)  an 
axisymmetric  shear  layer  originating  from  the  sharp  corner  of  the  base.  Tur¬ 
bulent,  subsonic  (incompressible)  shear  layers  are  characteristic  of  developing 
large  coherent  structures  that  are  of  a  highly  unsteady  nature.  (For  two  di¬ 
mensional  shear  layers,  these  large  structures  are  strongly  two-dimensionally 
organized.)  Thus,  it  is  not  surprising  that  large  coherent  structures,  which 
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can  be  observed  to  originate  in  the  shear  layers  directly  behind  the  blunt  base, 
are  present  also  for  subsonic  turbulent  wakes  behind  blunt  trailing  edges. 

There  is  considerable  evidence  that  the  cause  of  the  large  structures  is 
due  to  the  hydrodynamic  instability  of  the  (time-averaged)  mean  flow  and 
that  the  development  of  these  structures  can  be  captured  by  stability  the¬ 
ory.  In  fact,  certain  aspects  of  the  development  can  be  captured  even  with 
linear  stability  theory  although  intensities  (amplitudes)  of  the  structures  are 
often  too  large  for  the  linear  stability  theory  to  be  valid.  Experimental  re¬ 
sults  for  incompressible  turbulent  mixing  layers,  two-dimensional  turbulent 
wakes,  and  axisymmetric  wakes  with  a  blunt  base  and  comparison  with  lin¬ 
ear  stability  theory  have  shown  that  certain  key  features,  such  as  dominant 
frequencies,  mode  shapes  (amplitude  distributions),  and  streamwise  spacing 
(streamwise  wave  lengths)  of  the  structures  can  be  well  predicted  by  linear 
stability  theory  (Wygnanski  et  al.  (1986),  Marasli  et  al.  (1989)).  These  in¬ 
vestigations  support  the  notion  that  hydrodynamic  instabilities  give  rise  to 
the  generation  and  development  of  large  structures. 

The  dynamical  behavior  of  these  structures  is  responsible  for  the  strong 
unsteady  flow  behavior  in  the  wake.  Thus,  in  reality,  there  is  no  steady  tur¬ 
bulent  wake  flow  (and,  in  particular,  there  is  no  dead  wake  region  near  the 
base  of  a  bluff  body),  even  when  the  small  scale  (high  frequency)  fluctua¬ 
tions  are  not  taken  into  consideration.  Rather,  the  flow  is  highly  unsteady. 
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dominated  by  large  amplitude  and,  relative  to  the  small  scales,  low  frequency 
fluctuations.  The  steady  mean  flow  measured  in  experiments  is  in  fact  the 
time  average  of  the  time  dependent  fluctuating  flow  fleld.  Since  the  ampli¬ 
tude  of  these  large  scale  fluctuations  can  be  relatively  large,  the  total  flow 
field  cannot  actually  be  composed  by  linear  superposition  of  the  various  fluc¬ 
tuating  components  onto  the  (steady)  mean  flow  field.  Rather,  because  of 
nonlinear  interaction  between  the  various  fluctuation  components,  the  actual 
mean  flow  may  be  strongly  dependent  on  the  composition  of  the  fluctuating 
parts  of  the  flow  field. 

This  evidence  was  confirmed  experimentally  for  incompressible  wake  flows 
Wygnanski  et  al.  (1986),  Marasli  et  al.  (1989)  by  artificially  forcing  the  flow. 
With  artificial  forcing,  existing  modes  (structures)  could  be  modified  or  new 
modes  (structures)  created.  As  a  consequence,  the  mean  flow  could  be  mod¬ 
ified  substantially.  This  study  showed  convincingly  that  there  is  no  unique 
turbulent  mean  flow  for  the  wake  and  not  even  for  the  far  wake.  Rather,  the 
mean  (steady)  wake  flow  is  strongly  dependent  on  the  nature  of  the  large- 
scale  coherent  structures  and,  as  a  consequence,  on  their  dynamical  behavior. 

With  this  in  mind,  it  is  not  surprising  that  mean  flow  fields  depend 
strongly  on  experimental  conditions  and  can  even  vary  for  the  same  facil¬ 
ity  with  minor  changes  of  the  experimental  parameters.  Also,  for  this  rea¬ 
son,  it  is  not  surprising  that  current  turbulence  models  for  calculating  mean 


flows  using  the  Reynolds  averaged  Navier  Stokes  formulations  are  perform¬ 
ing  poorly  for  flows  that  are  known  to  be  dominated  by  the  dynamics  of 
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large  structures,  that  is,  for  flows  with  massive  separations  and  free  shear 
layers.  Chances  to  arrive  at  better  turbulence  models  are  rather  slim  un¬ 
less  the  physical  and  dynamical  behavior  of  the  large  structures  are  better 
understood  and  until  this  knowledge  can  be  implemented  in  the  turbulence 
modeling. 

Thus,  in  the  context  of  the  proposed  research,  the  question  arises;  Do 
large  structures  play  a  similarly  important  role  for  supersonic  separated  flows 
and  in  particular  for  supersonic  axisymmetric  wakes?  Unfortunately,  there 
are  relatively  few  experimental  investigations  that  have  focused  on  this  issue. 
However,  when  looking  at  flow  visualization  pictures  of  supersonic  wake  flows, 
very  distinct  patterns  with  large  scale  structures  can  be  observed. 

For  planar  supersonic  shear  layers,  the  experimental  investigations  of  Pa- 
pamoschou  &  Roshko  (1988)  have  convincingly  demonstrated  the  presence  of 
dominant  large  scale  structures  for  a  variety  of  conditions,  e.g.,  Mach  num¬ 
bers  and  density  ratios.  In  fact,  relative  to  the  shear  layer  thickness,  these 
structures  are  often  considerably  larger  (and  therefore  appear  to  be  more 
dominant)  than  those  for  incompressible  shear  layers.  Large  scale  (coher¬ 
ent)  structures  for  supersonic  shear  layers  were  also  observed  by  Ortwerth  & 
Shine  (1977)  and  for  supersonic  jets  by  Oertel  (1979).  For  the  supersonic  ax- 
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isymmetric  shear  layer  behind  a  backward-facing  step  investigated  by  Roshko 
&  Thomke  (1966),  flow  visualization  suggests  a  mean  flow  that  is  predomi¬ 
nantly  periodic  in  the  azimuthal  direction,  indicating  possibly  the  presence 
of  large  (coherent)  structures  with  a  dominant  azimuthal  wave  length  [see 
Figure  6  in  Roshko  &  Thomke  (1966)]. 

The  fact  that  wind  tunnel  interference  and  interference  from  model  sup¬ 
port  are  strongly  affecting  the  mean  flow  behavior  for  supersonic  axisymmet- 
ric  wake  flow  experiments  may  be  an  indication  that  this  might  be  caused 
by  the  presence  of  large  coherent  structures.  The  dynamical  behavior  of  the 
large  structures  may  be  changed  by  the  interference  disturbances,  or  addi¬ 
tional  large  structures  may  be  generated  due  to  changes  in  the  hydrodynamic 
stability  behavior  of  the  mean  flow,  as  discussed  above. 

Some  quantitative  evidence  of  the  existence  of  dominant  large  structures 
in  supersonic  axisymmetric  wake  flows  was  provided  by  the  experiments  of 
Demetriades  (1968),  who  investigated  the  unsteady  nature  of  the  flow  fleld. 
The  amplitude  spectra  (see  Figures  11  to  14  in  Demetriades)  display  distinct 
peaks  at  certain  (relatively  low)  frequencies,  thus  indicating  the  presence  of 
dominant  modes  (structures).  This  is  supported  by  plots  of  rms-amplitude 
distributions  (with  respect  to  radial  distance  from  the  centerline)  for  the  axial 
velocity  fluctuation  and  the  density  fluctuation  (Figures  2  and  3,  respectively, 
of  Demetriades).  The  amplitude  distribution  for  the  velocity  fluctuation  is 
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highly  reminiscent  of  that  of  an  incompressible  axisymmetric  wake,  where 
it  is  known  that  this  profile  is  due  to  the  presence  of  dominant,  large  co¬ 
herent  structures  (Cannon  &  Champagne  (1991)).  In  fact,  Morkovin  (1968) 
suggested  long  ago  that  when  normalized  properly,  the  distribution  of  fluctu¬ 
ations  (as  caused  by  large  structures)  has  a  universal  character,  even  in  very 
diverse  flow  regimes,  e.g.,  even  when  comparing  subsonic  and  supersonic 
flows. 

In  a  well-planned  and  carefully  executed  experimental  program  at  the 
University  of  Illinois,  Dutton  and  co-workers  [see,  e.g..  Smith  &:  Dutton 
(1996)]  identified  coherent  structures  for  a  supersonic  two-dimensional  wake 
behind  a  blunt  base.  The  observed  dominant  structures  appear  to  be  sig¬ 
nificant  relative  to  the  shear  layer  thicknesses.  These  observations  were  in 
qualitative  agreement  with  those  of  Papamoschou  &  Roshko  (1988)  for  a 
planar  supersonic  shear  layer.  From  our  previous  numerical  investigations 
of  transitional  and  turbulent  base  flows  (two-dimensional  and  axisymmetric; 
albeit  at  lower  Reynolds  numbers  than  the  experiments,  see  section  5),  we 
found  highly  energetic  structures  with  amplitudes  of  about  15%  of  the  free 
stream  velocity  that  may  have  a  significant  effect  on  the  global  flow  behavior 
(see  earlier  discussion). 

More  recently,  employing  sophisticated  flow  visualization  techniques,  Dut¬ 
ton  and  co-workers  (Bourdon  &;  Dutton  (1998,  20006))  were  able  to  clearly 


show  the  presence  of  large  turbulence  structures  in  the  axisymmetric  base 
flow.  Therefore,  from  our  previous  simulations  and  the  recent  experiments 
by  Dutton  and  co-workers,  we  have  sufficient  evidence  that  large  coherent 
structures  are  present  in  supersonic  axisymmetric  base  flows  and  may  in¬ 
deed  play  an  important  role.  The  fact  that  dominant  structures  can  exist  in 
axisymmetric  wake  flows  and  may  play  an  important  role  is  crucial  to  the 
proposed  research,  as  numerical  simulations  will  be  employed  for  these  inves¬ 
tigations.  For  numerical  simulations  to  be  useful  for  providing  new  insight 
into  the  physics  of  supersonic  base  flows,  it  is  therefore  essential  that  the 
simulation  techniques  are  such  that  the  role  of  the  large  turbulent  structures 
is  accurately  represented  and  the  dynamics  of  the  structures  can  be  captured 


reliably. 


32 


3.  GOVERNING  EQUATIONS 

The  flow  is  assumed  to  be  an  ideal  gas  with  constant  specific  heat  coeffi¬ 
cients  and  the  following  assumptions  are  made  for  the  constitutive  relations: 
the  flow  can  be  treated  as  a  Newtonian  fluid  with  zero  bulk  viscosity,  the 
Fourier  law  is  valid  for  heat  conduction,  the  Prandtl  number  is  constant  in 
the  region  of  interest  and  the  viscosity  can  be  evaluated  according  to  Suther¬ 
land’s  law. 

3.1  Dimensional  Equations 

The  governing  field  equations  for  this  research  are  the  three-dimensional 
compressible  Navier-Stokes  equations  in  cylindrical  coordinates.  Body  forces 
are  neglected  and  the  Prandtl  number  is  assumed  constant.  The  conservative 
form  of  the  equations  was  chosen  because  it  can  be  put  in  vector  form.  This 
method  was  used  by  Thumm  (1991)  and  is  described  in  Anderson  et  al. 
(1984).  Following  are  the  equations  in  tensor  notation,  with  an  asterisk 
denoting  a  dimensional  variable. 

Continuity: 

Pr  +  =  0  (3.1) 

Momentum: 


ip-u’)  ,.  +  (p-u‘ul  +p'5,t  -  4),^.  =  0 


(3.2) 


Total  energy: 


(p-E-),,.  +  {p-ulH-  +  gj  -  =  0  (3.3) 

where  the  total  energy  and  the  total  enthalpy  are 

E*  =  c*r*  +  rWjuj  and  if*  =  £■*  +  ^  ,  respectively.  (3.4) 

/  p* 

The  stress  tensor  is 


'^ik  2/i  //  with  —  -  {u*k*  (3.5) 


and  the  heat  flux  vector 

ql  =  -k*T%.  .  (3.6) 

Finally,  the  equation  of  state  closes  the  system  and  is  given  as 


p*  =  p*TVT* . 


Sutherland’s  Law  is  used  to  calculate  the  viscosity: 


o* 

I 

Uv  ' 

KT*  +  n*sJ 

with 


TVsu  =  iio.eif 

3.2  Non-dimensional  Equations 


(3.7) 


(3.8) 


(3.9) 


All  quantities  are  made  dimensionless  by  conditions  at  a  reference  point 
in  the  flow,  which  in  this  case  is  the  freestream/inflow  point,  and  a  reference 
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length,  which  will  be  chosen  to  be  the  body  radius.  All  field  equations  and 
constitutive  relations  are  made  dimensionless  via  the  set  T^,  r*,  k*^. 

This  leads  to  the  following  dimensionless  groups: 


Re=&l^,Pr-!^ 


Ec= 


(3.10) 


^p-^OO 


Mr^  = 


% 

a* 


The  mass,  momentum  and  energy  equations  are  made  dimensionless  by  mul¬ 
tiplying  by  r*/p*^u*^,  r*/p^w^,  and  r*/p*^u^,  respectively.  This  leads  to 
the  non-dimensional  equations,  now  given  in  vector  form  in  cylindrical  coor¬ 
dinates. 


with 


m  dA 
dt  ^  dz  ^ 


dB  IdC  1^  „ 

~  ^ 

or  r  ou  r 


(  P\ 


U  = 


pu 

pv 


pw 

\pEj 


( 


A  = 


pu  \ 

+  p  -  r22 
PUV  -  Trz 


PUW  -  Tgz  I 

\  puH  +  qz-  UTzz  -  VTrz  -  WTgz  ) 


B  = 


pv 

pUV  -  Trz 
pv"^  +  p-  Trr 
PVW  -  TrO 

pvH  +  Qr  —  UTrz  —  VTrr 


\ 


WTr,  } 


(3.11) 


(3.12) 


(3.13) 


(3.14) 
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/ 


pw 


\ 


C  = 


D  = 


puw  -  Tez 
PVW  -  Tr  theta 
pw^  +p- Tee 

pwH  +  -  UTez  -  vTre  -  wree  J 

/  /w  \ 

puv  -  Trz 

pv^  -  pw^  -  Trr  +  Tee 

2pvw  —  2Tre 

pvH  +  qr-  UTrz  “  VTrr  -  WTre  J 


with  the  stress  components 


Tzz 

Trr 

Tee 

'^rz 


Tez 


Tre 


1 

Re 

1 

Re 

_1_ 

Re 


2p,Szz 

2p,Srr 

2p,See 


Re 


2pSre 


-p  {Szz  +  Srr  +  See) 
2 

gM  {Szz  +  Srr  +  See) 
2 

{Szz  +  Srr  +  See) 


The  strainrates  in  cylindrical  coordinates  are: 


_  du  n  _  c  _  ^  ^ 

~  >  Srr  —  )  See  —  H 

oz  or  r  oU  r 


(3.15) 


(3.16) 


(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 


(3.23) 


1 

2 


/  du 

V  ^  dz) 


and 


1  (Idv  d  w\ 
2\rd9^^drr) 

(3.24) 
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Non-dimensionalizing  Sutherland’s  Law  yields 

Note  that  since  the  Prandtl  number  is  assumed  constant  throughout  the  flow 
field,  the  dimensionless  thermal  conductivity  and  viscosity  are  identical: 


3.3  Averaged  Equations 


Following  Speziale  (1997),  applying  a  filter  to  the  governing  equations, 
the  fundamental  equations  for  LES  are  found: 
the  filtered  continuity  equation 

the  filtered  momentum  equations 

d  d 

M  ^  'dx'k  “  p(^ik)]  =  0  ,  (3.28) 


and  the  resolved  energy  equation 


—  (pEr)  +  —  [pUkHR  +  qk  +  Qk-Ui  (Tik  -  paik)]  =  11 , 


(3.29) 


where  the  overbar  represents  a  standard  filter  (i.e.,  for  this  research,  the 
implicit  filter  inherent  to  the  discretization  scheme  with  a  filter  width  related 
to  the  grid  size)  and  a  tilde  denotes  the  Favre-average  defined  as 
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7=^.  (3-30) 

P 

The  velocity  vector  Ui  and  the  temperature  T  are  decomposed  into  a  Favre- 
average  and  a  fluctuating  part  (denoted  by  ”)  and  the  pressure  p  and  the 
density  p  are  split  into  the  Reynolds-average  (denoted  by  -)  and  a  fluctuating 
part  (denoted  by  ’). 

The  resolved  energy  E,  the  resolved  total  enthalpy  the  resolved  heat 
flux  vector  qk  and  the  resolved  stress  tensor  Tik  are 

E  =  CyT  4-  —UjUj  ,  H  =  E  +  —  ,  Qk  —  ^ PrEcRe ) 

Tik  =  (^ik  -  ^SjjSii^  with  Sik  =  -  {ui,k  +  Uk,i)  ,  (3.32) 

respectively.  The  resolved  molecular  viscosity  is  computed  according  to 
Sutherland’s  law  (see  e.g.  White  (1991)) 


(t)  = -c  (f )  =  f 


with  = 

T  +  TlsuJ 


(3.33) 


Finally,  the  non-dimensional  equation  of  state  determines  the  resolved 


pressure 


p  = 


W 


(3.34) 


For  the  current  work,  7  =  1.4  and  Pr  =  0.70. 
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The  above  equations  contain  three  terms  that  do  not  occur  in  the  unfil¬ 
tered  equations:  the  subgrid  stress-tensor  aik,  the  subgrid  heat-flux  vector 
Qk  and  the  source  term  11  in  the  energy  equation.  These  terms  have  to 
modelled  (see  section  3.4).  In  case  of  a  DNS,  where  one  assumes  that  all 
relevant  time  and  length  scales  are  resolved  by  the  computational  grid,  the 
model  terms  are  set  to  zero  implying  <j)  —  <f>  =  In  the  other  limit,  where 
the  filter-width  is  so  large  that  all  fluctuations  are  filtered  out,  a  traditional 
RANS  is  recovered. 

3.4  Turbulence  Models 


If  the  numerical  simulations  are  not  performed  in  the  DNS  limit,  a  closure 
is  needed  for  equations  (3.27)  -  (3.29). 

The  subgrid  heat-flux  vector  is  modelled  by  assuming  similarity  in  the 
gradient  transport  of  heat  and  momentum,  relating  the  thermal  eddy  dif- 
fusivity  kt  to  the  eddy  diffusivity  nr  (see  below)  by  a  constant  turbulent 
Prandtl  number  Ptt  as  described  in  Speziale  &  So  (1996) 


Qk 


'yEc  Pr-r 


Tk. 


(3.35) 


Speziale  &  So  (1996)  suggest  that  for  most  engineering  applications  a  tur¬ 
bulent  Prandtl  number  of  Pry  =  0.9  is  adequate.  Note  that  for  a  turbulent 
Prandtl  number  of  unity  complete  similarity  is  assumed  which  is  the  well 


known  Reynolds  analogy  for  turbulent  heat  transfer. 
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The  additional  source  term  11  in  the  energy  equation  includes  a  pressure 
dilatation  term,  and  terms  involving  turbulent  dissipation  and  the  subgrid 
mass-flux.  It  was  shown  by  Sarkar  et  al.  (1991)  that  not  only  the  pres¬ 
sure  dilatation  term  but  also  the  compressible  dissipation  are  important  in 
compressible  turbulence.  Therefore,  both  effects  are  modelled  according  to 
Sarkar  et  al.  (1991)  and  Sarkar  (1992) 

f/u'-  j  =  a^paikSikMT  +  a^peM^ ,  (3.36) 


with  the  turbulent  Mach  number  (a  measure  of  influence  of  compressibility 
on  turbulence) 


2K* 


(3.37) 


The  turbulent  dissipation  is  defined  as  and  the  subgrid  mass-flux 

is  modelled  as  u'-  =  -r—p-Tpi  •  Combining  these  terms  yields  the  extra  source 
p(^P 


term 


,k 

(3.38) 

with  the  constants  02  =  0.15  and  03  =  0.2. 

For  the  computation  of  the  subgrid-stress-tensor  atk,  an  explicit  alge¬ 
braic  stress  model  is  used.  It  was  derived  from  the  explicit  solution  of  the 
equilibrium  form  of  the  modelled  Reynolds  stress  transport  equation  for  in¬ 
compressible  flows  in  Gatski  &  Speziale  (1993)  and  extended  to  compressible 


n  =  (1  -  02 Mt)  pcTikSik  +  (1  -  asMl)  pe  -  (r^fe  -  p6ik) 


/xr  /I 


flows  by  Speziale  (1997).  The  model  is  of  the  form  of  an  anisotropic  eddy 
viscosity  model  with  strain-dependent  coefficients: 
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=  \KSik  -  (Si»  - 

+  {SijSjk  -  -SpgSpgSikj  .  (3.39) 


To  remove  the  possibility  of  singularities,  a  regularized  version  for  a*  is  used: 


3 

-  2rp  6^2 


3(1+7?^) 

*3  -H  772  +  6772^2  +  6^2  ■ 


(3.40) 


The  terms  77  and  ^  are  related  to  the  irrotational  and  the  rotational  strain 
rate  invariants  according  to 


"  2S77 


(3.41) 

(3.42) 


For  compressible  calculations  the  constants  ccj  are  given  as  ai  =  0.37436, 
0(2  =  0.11518,  and  013  =  0.10799.  For  ’’incompressible”  boundary  layer  test- 
calculations  the  incompressible  constants  a\  =  0.227,  a2  =  0.0423  and  0:3  = 
0.0396  were  used^ .  The  turbulent  eddy  viscosity  is  given  by 


pK^ 

fJ'T  =  with  Cp  =  0.09  (3.43) 

The  turbulent  kinetic  energy  K  and  the  turbulent  dissipation  rate  e  are 
computed  solving  two  additional  transport  equations  (see  Wilcox  (1998)) 


^For  Qi  =  2Cp  and  0:2  =  03  =  0  equation  (3.39)  reduces  to  the  standard  K  -e  model. 
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which  are  given  here  in  dimensionless  form: 


_  ^  4  _ 

+  CezpRe^—  -  -peui,i ,  (3.45) 

with  the  constants  ck  =  1-0,  =  1.3,  C^i  —  1.44,  0^2  =  1-83  and  C^z  = 

0.001.  The  source  term  of  eqn.  3.44  is  the  negative  value  of  the  source  term 
in  the  energy  equation  because  any  production  of  turbulence  kinetic  energy 
must  be  equal  to  the  dissipation  of  the  averaged  energy.  The  turbulent 
Reynolds  number  is  defined  as: 


Rer  =  ^  .  (3.46) 

pe 

To  remove  a  singularity  at  walls  (i.e.,  R!"  =  0)  in  the  destruction  term  of  the 
^-equation,  a  damping  function  of  the  form 


/e2  =  1  —  exp  ^—ReVO.lKN^ 


(3.47) 


is  used,  where  N  is  the  wall-normal  distance.  In  the  ASM  model,  this  wall¬ 
damping  function  is  the  only  term  containing  a  wall  distance,  the  Reynolds 
stress  model  automatically  accounts  for  near  wall  effects  through  the  compu¬ 
tation  of  Oi\.  To  be  completely  independent  of  any  kind  of  wall  function,  the 
idea  of  Durbin  (1993)  is  also  implemented.  It  incorporates  the  assumption 


that  the  smallest  physical  time-scale  in  a  turbulent  flow  is  the  Kolmogorov 
time-scale.  Therefore,  by  computing  the  turbulence  time-scale  as 
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(3.48) 


with  Cx  «  6  and  setting  /£2  =  1,  a  completely  wall-distance  independent 
model  is  obtained. 


3.5  Flow  Simulation  Methodology  (FSM) 


An  important  aspect  of  ongoing  research  is  the  application  of  the  Flow 
Simulation  Methodology  (FSM).  The  fundamental  concept  of  the  methodol¬ 
ogy  was  already  introduced  in  section  1.  The  details  of  the  implementation 
are  given  now. 

In  FSM  the  turbulent  stress  tensor  that  is  fed  back  to  the  filtered  Navier- 
Stokes  equation  is  scaled  with  a  contribution  function, 


cri,  =  /(A/L,)<Tj,  (3.49) 

where  A  =  [{Ax‘^  +  Ay^  +  A2:^)/3]^/^  is  the  representative  computational 
grid  size  and  Lk  =  is  the  Kolmogorov  length  scale.  The  final  form 

of  /(A/Lfc)  has  not  yet  been  determined,  but  for  the  current  work  it  was 
chosen  to  be  of  the  form  proposed  by  Speziale  (1997): 


( 

/(A/Lfc)  =  1  -  e  Lk 


n 


(3.50) 
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where  /?  and  n  are  adjustable  parameters. 

As  the  ratio  ^  becomes  small  (the  grid  resolution  is  sufficient  to  re¬ 
solve  structures  the  size  of  the  Kolmogorov  scale),  /(A/L*)  approaches  zero 
and  the  computation  will  approach  the  DNS-limit.  For  insufficient  resolu¬ 
tion,  /(A/Lfc)  approaches  unity  and  the  RANS-limit  is  approached.  For  all 
intermediate  values  of  the  contribution  function,  a  ’’non-traditional”  Large 
Eddy  Simulation  (LES)  based  on  a  state-of-the-art  Reynolds  Stress  model  is 
performed. 
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4.  NUMERICAL  PROCEDURE 
4.1  Overview 

For  the  spatial  derivatives,  we  are  now  using  high  order-accurate  finite- 
difference  approximations  (sixth  order  in  the  axial  and  radial  directions)  and 
a  pseudo-spectral  (Fourier)  decomposition  in  the  azimuthal  direction.  This 
high  accuracy  is  necessary  for  efficiently  carrying  out  the  highly  demanding 
simulations  of  the  dynamical  behavior  of  the  large  structures  in  supersonic 
wake  flows.  In  fact,  during  the  years  with  funding  from  the  ARO  grant,  we 
have  drastically  improved  the  spatial  accuracy  of  our  Navier-Stokes  code. 
Direct  Numerical  Simulations  of  turbulent  wake  flows  for  higher  Reynolds 
numbers,  as  performed  in  this  research,  require  increasingly  better  numerical 
resolution  for  increasing  Reynolds  numbers  as  the  size  of  the  spatial  scales 
decreases  with  Reynolds  number.  The  use  of  finer  and  finer  grids  is  limited 
by  the  memory  and  computing  power  of  currently  available  supercomputers. 

Therefore,  using  spatial  difference  approximations  with  higher  accuracy 
will  allow  numerical  resolution  of  the  increasingly  smaller  scales  without  re¬ 
quiring  a  decrease  of  the  grid  intervals.  From  experience  with  our  previous 
fourth-order-accurate  code  (using  standard  difference  approximations)  and 
from  our  experience  using  compact  finite  differences  for  incompressible  sim¬ 
ulations  (Meitz  &  Fasel  (2000)),  we  came  to  the  conclusion  that  the  consid- 
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erable  demands  regarding  accuracy  for  simulations  of  turbulent  supersonic 
wakes  could  be  best  met  by  employing  compact  difference  approximations  of 
at  least  fourth-order  accuracy.  Toward  this  end,  we  have  entirely  revamped 
our  previously  developed  code  and  implemented  compact  differences  where 
the  order  of  accuracy  is  adjustable  between  fourth  and  sixth  order.  This  is  ac¬ 
complished  by  allowing  the  coefficients  in  the  approximations  to  be  adjusted 
so  that  the  accuracy  (spectral  accuracy)  can  be  tailored  to  the  particular 
needs  of  the  simulation. 

The  advantage  of  employing  high-order  compact  differences  versus  stan¬ 
dard  differences  is  apparent  from  Figure  4.1,  where  the  numerically  modified 
wave  number  of  the  computed  scales  is  compared  with  the  actual  (physical) 
wave  number  of  the  scales  (which  should  be  resolved  in  the  simulation).  It  is 
obvious  that  the  drop-off  in  these  modified  wave  numbers  starts  much  earlier 
for  a  standard  fourth-order  approximation  than  for  a  compact  difference  of 
the  same  formal  accuracy. 

Use  of  methods  with  high  accuracy,  and  thus  employing  compact  differ¬ 
ences,  is  equally  beneficial  for  Large-Eddy  Simulations  (LES)  and,  in  partic¬ 
ular,  for  our  new  FSM  (see  discussion  in  3.5).  Of  course,  the  higher  accuracy 
allows  use  of  larger  step  sizes  for  the  Navier-Stokes  computations  for  the 
resolved  scales  and  therefore  reduces  demands  on  computer  memory  and 
computing  power.  In  addition,  and  this  is  crucial,  the  higher  accuracy  allows 
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Wavenumber 

Figure  4.1  Modified  wave  number  versus  physical  wave  number.  Comparison 
of  standard  and  compact  difference  approximations. 

rigorous  testing  and  validation  of  the  turbulence  models  (required  for  the 
unresolved  subgrid  scales),  which  have  to  be  employed  in  LES.  Using  lower- 
order-accurate  difference  approximations  does  not  allow  a  clean  separation 
of  the  effects  of  the  subgrid-scale  turbulence  model  from  the  eflFects  of  the 
truncation  error  and,  therefore,  the  role  of  the  model  cannot  be  scrutinized. 
In  fact,  in  many  LES  investigations  published  in  the  literature,  the  trunca¬ 
tion  error  of  the  numerical  method  played  a  more  important  role  than  the 
subgrid-scale  model  itself  and  therefore  totally  contaminated  the  results. 
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4.2  Vectorized  Equations  in  Cylindrical  Coordinates 


All  transport  equations  given  in  section  3  can  be  written  in  the  form  of 
an  advection  and  a  diffusion  term  on  the  left  hand  side  and  a  source  term  on 
the  right  hand  side.  In  cylindrical  coordinates  the  structure  is: 


dU  dA  dB  IdC  1^  ^ 

- h  - h  +  -D  —  S  . 

at  oz  or  r  ad  r 
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xStJ 

with  n,  XI K  and  He  being  the  source  terms  given  in  equations  3.38,  3.44  and 
3.45  respectively. 


4.3  Azimuthal  Fourier  Transforms 


Water  tunnel  results  conducted  by  Siegel  (1999)  show  that  the  flow  picks  a 
symmetry  plane  that  moves  on  a  time-scale  much  slower  than  that  of  interest 
for  the  unsteady  structures  in  subsonic,  incompressible  flows.  As  in  previous 
work  on  this  topic  by  Tourbier  (1996),  symmetric  Fourier  transforms  are 
therefore  used  in  the  azimuthal  direction  to  achieve  a  signiflcant  reduction  in 


49 


computing  time.  Previous  work  by  Bourdon  &  Dutton  (1998)  shows  that  even 
the  compressible  flow  shows  some  symmetry  in  azimuthal  direction  and  the 
most  recent  publication  by  Bourdon  &  Dutton  (2001)  shows  a  clear  symmetry 
of  the  flow  in  circumferential  direction  if  delta-shaped  disturbances  on  the 
body  are  used.  Note  that  using  symmetric  transforms  enforces  a  symmetry 
plane  similar  to  the  way  the  surface  disturbances  do  in  the  experiment.  For 
simplicity,  only  the  equations  for  the  DNS  version  of  code  (5  equations  solved 
in  vector  form)  are  shown  in  Fourier  space. 

Replacing  each  flow  variable  with  the  appropriate  sine-  or  cosine-  trans¬ 


form  gives 
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Y' 
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A'l{r,  z)  sin  k9 
yAl{r,  z)  cos  k9  J 

(  C{(r,  z)  cos  k9 
£)^(r,  2)  cos  k9 
b\{r,  z)  cos  k9 
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k=0 
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/Bf(  r,  z)  cos  k9\ 
B^ir,  z)  cos  k9 
B^{r,  z)  cos  k9 
B^{r,  z)  sin  k9 
\B^{r,  z)  cos  k9  J 


(4.8) 


Commuting  derivative  and  summation,  and  using  the  orthogonality  of  the 


sine  and  cosine  function  yields  an  equation  for  each  Fourier  mode  k, 
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/  kC^{r,z)  +  D![{r,z)  \ 
kCj{r,z)  +  Dlir,z) 

H —  kCt(r,z)  +  D^{r,z)  =0, 

-kCl{r,z)  +  Dl{r,z) 

\  kCl{r,z)  +  DUr,z) 

(4.9) 

k  =  0,...,K. 

This  set  of  equations  is  now  solved  with  finite  differences  in  the  radial  and 
streamwise  direction  as  described  in  section  4.5. 

4.4  Parity  Conditions  for  Fourier  Modes  of  Functions  in  Cylin¬ 
drical  Coordinates 


UHr,z) 


z) 


B^{r,  z) 


g  Ujir,z)  Q  A}{t,z)  q  Bj{r,z) 


UHr,z)  \+^  A\{r,z)  + 
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U^{r,z) 


Ai(r,z) 
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Bjir,  z) 
Blir,  z) 
Blir,  z) 


Another  challenge  of  this  research  has  been  the  computation  at  the  axis, 
which  in  cylindrical  coordinates  is  represented  by  r  =  0.  The  treatment  of 
this  boundary  is  discussed  in  the  following.  Axis  boundary  conditions  are 
derived  for  the  Fourier  modes  of  functions  in  cylindrical  coordinates  that  are 
necessary  for  the  function  to  be  smooth. 

Consider  smooth  (C°°)  functions,  either  scalar,  vector  or  second  order 
tensor  in  a  cartesian  coordinates  s,  t  rotated  by  angle  (j).  Smooth  implies  a 
convergent  Fourier  representation, 

OO 

/(>■.«)  = 

fe=0 

The  complex  representation  is  used  only  for  convenience,  as  all  functions 
of  interest  here  are  real- valued. 
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On  s,  the  cylindrical  coordinates  and  unit  vectors  are  given  by, 


r  =  |s 


(4.11) 


s^o-  ds'^  i^o+  ds” 


(4.15) 


Even  if  /  is  a  smooth  function,  radial  derivatives  of  f  need  not  be  con¬ 
tinuous  at  the  origin  due  to  the  coordinate  singularity.  Even  so,  the  Fourier 


modes  of  functions  that  are  the  radial  derivatives  of  smooth  functions  have 


notable  parity  properties.  These  are  seen  almost  by  inspection  of  the  Fourier 
representation,  given  that  fkij)  is  either  even  or  odd  in  r.  One  may  note  that 
any  odd  order  of  radial  differentation  will  switch  the  parity  ’’state”  relative 
to  the  original  function,  while  an  even  number  will  not. 

Define  ’’even  parity”  as  being  an  even  function  in  r  in  the  even  Fourier 
modes,  while  being  an  odd  function  in  r  in  the  odd  modes.  ’’Odd  parity”  is 
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then  being  an  even  function  in  the  odd  modes,  etc.  For  functions  with  even 
parity, 


even  k-\-n  even 
odd  k-\-n  odd 


Functions  with  odd  parity  have 


even 

odd 


k-¥n  odd 
k  -\-n  even 


Note  that  derivatives  in  6  do  not  change  the  parity  of  a  function. 


(4.16) 


(4.17) 


4.4.1  Scalars 


The  Fourier  representation  of  /  is 

oo 

/(^)  =  A(kl)e*^^sgn''(s) ,  (4.18) 

*=0 

or  in  general. 


=  5^/fe"^(|s|)e'^‘^sgn''(s)sgn"(s) .  (4.19) 

k=Q 

Substituting  this  into  equation  4.15, 

OO  OO 

l]/fe”^(0)e**^(-l)(''+”)  =  (4.20) 

fc=0  k=0 


00 

Y  [1  -  (-1)^'=+")]  =  0  (4.21) 

A:=0 

oo 

Y  2/f)(0)e*'=^  =  0. 

k—OyTi-i-k  odd 


(4.22) 


Since  this  is  required  for  all  <f>, 


/i“>(0)  =  0,Vn  +  *:odd, 


(4.23) 


that  is,  scalar  functions  have  even  parity. 

4.4.2  Vectors 


Consider  now  a  smooth  vector  field  f  with  (convergent)  Fourier  represen¬ 
tation 


(X) 

f  (r,  2)  =  ^/fc(r,  z)er  +  gk{r,  z)ee  -f  hk{r,  z)e^  e**®  (4.24) 


fc=0 


with 


e.  =  sgn(s)e. 
=  sgii(s)e, 


(4.25) 

(4.26) 


On  s,  f  is  represented  by, 


=  [(A(l«l>-2^)es+5fe(|5U)e«)sgn(s)  +  hfc(|s|,2)e^  e*''^sgn*(s) 

A;=0 

(4.27) 

Taking  the  n-th  derivative  in  s  and  taking  the  two-sided  limit  as  before, 
we  obtain 


00 

»=i: 


< 


fc=0 
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+^ifc,nr(0,2)ef  [1  -  (-1)"+*=+!] 
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(4.28) 
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or, 

OO 

0=  2A,„(0,  z)e.e‘W  (4.29) 

n+fe  even 

OO 

+  2^A:,nr(0, 2)646**''^ 

k=Q 

71+ k  even 

OO 

+  Y1  2AMr(0, -2)6.6'*=^ 

fc=0 

n+k  odd 

Again  invoking  orthogonality, 

/fc,nr(0,  z)=0,n-\-k  even 

^fc,nr(0,  z)  =  0,  n  +  k  even  (4.30) 

hk,nr{0,  z)  =0,  n  +  k  odd, 

so  the  radial  and  azimuthal  vector  components  have  odd  parity,  while  the 
axial  component  has  even  parity. 

4.5  Compact  Differencing  on  Arbitrary  Grids 

The  transport  equations  given  in  Fourier  space  in  equation  4.3  are  dis¬ 
cretized  with  compact  finite  differences  in  the  radial  and  streamwise  direc¬ 
tion.  Sixth-order  split  compact  finite  difference  stencils  are  derived  in  the 
following.  These  differences  are  valid  on  a  non-equidistant  grid  instead  of 
the  usual  mapping  of  the  computational  domain  to  a  uniform  grid.  This 
results  in  greater  flexibility  of  the  choice  of  grids  and  was  shown  to  lead  to  a 
lower  truncation  error  by  Meitz  &  Fasel  (2000). 
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4.5.1  Derivation 


Consider  an  ordered  set  of  points  {xj},  with  Xi  ^  Xj  for  i  ^  j.  It  is 
helpful  to  define  a  characteristic  grid  spacing  Ax  which  can  be  set  as  desired; 
Ax  =  xi—xo  is  chosen  here.  Note  that  this  implies,  mathematically  at  least, 
that  the  grid  is  refined  uniformly.  Define  a  dimensionless  grid  spacing  by 


4  = 


Xk  -  Xq 
Ax 


(4.31) 


Define  the  finite  difference  stencil  by 


^  ‘  Sx 

k=ki 


1 

t=tl 


(4.32) 


where  it  is  assumed  for  simplicity  that  the  stencil  spans  xq  for  both  function 
values  and  derivatives,  that  is. 


k\k2  <  0 


(4.33) 


I1I2  <  0 


The  stencil  is  homogeneous  in  its  coefficients.  One  condition  on  the  coeffi¬ 
cients  can  be  selected  arbitrarily. 

fca 

^  Ofc  =  1  (4.34) 

k=ki 

This  is  selected  to  allow  a  direct  comparison  of  the  truncation  error  term 
between  compact  and  standard  stencils.  The  expression  for  the  truncation 
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error  is 


=  («6) 
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Expand  f'{xk)  in  Taylor  series  about  xq-  Since  0°  is  not  defined,  first  pull 
out  the  zero  term  in  each  sum. 
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Exchange  order  of  summation  and  shift  an  index 
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Introduce  dimensionless  grid  and  collect  like  powers  of  Ax 
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Note  finally  that  Jq  =  0 
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A  diflference  formula  of  order  m  requires  all  terms  of  order  n  <  m  in  Ax 
to  be  zero.  Including  the  scaling  equation,  this  gives, 
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As  expected,  a  stencil  of  order  m  requires  m+2  equations.  This  is  the  same  as 
the  number  of  known  values  in  the  stencil,  (function  value  and  derivatives). 
When  these  conditions  are  met,  the  truncation  error  becomes, 


T  =  Arc" 
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(m+l) 


(m  +  1)! 


*2 


(m  +  l)afe5j 

.l=ll  k=^ki 


k=k\ 


(4.49) 
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This  can  be  arranged  in  a  matrix  equation  and  solved  for  specific  grid  point 
values.  Alternately,  the  matrix  can  be  inverted  symbolically  using  a  computer 
algebra  system  such  as  Maple  or  Mathematica,  and  the  solution  translated 
directly  into  source  code. 

4.5.2  Analysis  of  the  compact  diflFerence  stencils 

4.5.2. 1  Modified  Wavenumber 

Examine  the  effect  of  equation  4.32  when  applied  to  a  single  mode  f{x)  = 
giKx  exact  derivative  is 


(4.50) 


where  k  will  be  called  the  exact  wavenumber.  Define  the  modified  wavenum¬ 


ber  k'  by 


Sx 


(4.51) 


k'  is  found  by  substituting  4.51  into  4.32. 


k'Ax  =  —i 


"'^2  olKAxSk 


(4.52) 


4. 5. 2. 2  Amplitude  and  Phase  Error 


Consider  the  linear  model  convection  equation  in  one  dimension, 


(4.53) 
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with  initial  condition 


/(x,0)=e‘”* 

(4.54) 

The  exact  solution  at  t  =  At  is 

f{x,  At)  = 

(4.55) 

_  gZKXg““2i/(«Aa:) 

(4.56) 

where  u  is  the  CFL  number.  We  characterize  this  by  the  exact  amplitude 
and  phase,  given  by 


A  =  =  1 


(4.57) 


$  =  —ukAx  (4.58) 

If  the  standard  fourth  order  Runge-Kutta  scheme  is  applied  to  equation 
4.53  along  with  an  approximate  x-derivative, 

4  .. 

f{x,  At)  =  —  {—iuK'Ax)'^  .  (4.59) 

ik=ro  ■ 

The  numerically  obtained  amplitude  and  phase  are  then 


=  C  ^  (-w/c'Aa:)'' 


fe=0 


(4.60) 
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It! 


(4.61) 
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Figure  4.2  Radially  and  axially  stretched  grid  for  typical  DNS  (full  domain 
extends  to  r  =  10  and  z  =  25,  only  every  other  grid  point  shown). 


Examples  are  shown  for  an  exponentially  stretched  grid  and  various  com¬ 
pact  stencils  in  figures  4.3  and  4.4.  A  typical  grid  for  a  DNS  is  shown  in 
figure  4.2. 

Figures  4.3  and  4.4  show  an  analysis  of  different  spatial  discretizations 
in  combination  with  a  fourth-order  Runge-Kutta  method  at  CFL  =  0.15 
applied  to  the  linear  wave  equation.  It  becomes  clear  that  the  sixth-order 
compact  differences  have  a  far  superior  wave-number  accuracy  (see  Modified 
Wavenumber  plot).  However,  when  the  stretching  of  the  grid  is  accounted 
for  in  the  analysis,  the  amplitude  plot  shows  that  for  stretching  factors  larger 
than  unity,  amplification  will  occur  in  the  higher  wave-number  region.  Ther- 
fore,  a  spatial  filter  has  to  be  employed.  For  this  research,  a  fourth-order 


compact  filter  was  chosen. 
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Modified  Wavenumber  Stretching  Factor 


Amplitude  Normalized  Phase 


Figure  4.3  Modified  wavenumber,  Amplitude  and  Phase  distribution  for  a 
fourth  order  explicit  scheme  on  a  stretched  grid  in  combination  with  a  fourth- 
order  Runke-Kutta  time-integration  for  CFL  =  0.15 
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Figure  4.4  Modified  wavenumber,  Amplitude  and  Phase  distribution  for  a 
sixth  order  compact  scheme  on  a  stretched  grid  in  combination  with  a  fourth- 
order  Runke-Kutta  time-integration  for  CFL  =  0.15 
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4.6  Time  Integration 

A  fourth-order  Runge-Kutta  method  is  used  for  the  time-discretization. 
In  this  case,  four  intermediate  timesteps  are  used,  which  Ferziger  (1981) 
writes  as 


01  =  0n  +  -^/(0n) 

02  =  0n+— /(0i)  (4.62) 

03  =  0n  +  ^tf{(j>2) 

0n+l  =  0n  +  “^  [  /(  0n  )  +  2/(0i)  +  2/(^2)  +  /(03)]  • 


If  the  first  three  steps  are  substituted  into  the  last,  the  result  is 


=  ,!>„  +  At/M  +  ■  (4.63) 


Thus,  the  method  is  of  fourth  order  in  At.  As  written  above,  three  interme¬ 
diate  levels  of  variables  {4>i,  <f)2,  0$)  are  needed.  With  some  manipulation, 
the  form  used  in  the  code  can  be  obtained,  which  uses  only  two  intermediate 
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variables: 

,  ,  At  .  . 

Yl  —  0n  4  ^  Jyrn) 

At  ,  . 

<P2  =  (Pn  + 

4>i  =  +  2(f>2  (4.64) 

(t>2  =  <f>n  +  Atf{(f)2) 

^1  =  2  {—<i>n  +  01  +  02) 

0n+l  =  01  +  -^/(02)  • 

4.7  Integration  Domain 

For  the  numerical  simulations,  the  integration  domain  is  divided  into  two 
sub-domains.  One  domain  contains  the  boundary  layer  on  the  axisymmetric 
body  and  extends  to  the  corner  of  the  base.  The  second  domain  starts  at  the 
base  and  extends  to  the  outflow  boundary  (Figure  4.5).  As  in  the  simulations 
by  Tourbier  (1996)  only  the  last  part  of  the  approach  flow  is  computed  with 
the  inflow  boundary  layer  being  another  parameter  in  the  simulation. 

The  numerical  simulations  are  performed  in  two  steps.  In  the  first  step, 
a  steady  base  flow  is  computed  by  solving  the  Navier  Stokes  equation  for  an 
axisymmetric  geometry.  In  the  second  step,  the  unsteady,  three  dimensional 
flow  is  computed  using  the  steady  base  flow  (computed  in  step  1)  as  the  initial 
condition.  Thus,  the  Navier-Stokes  code  consists  of  two  main  components. 
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freestream  boundary 


Figure  4.5  Computational  Domain. 

one  for  calculating  the  steady,  axisymmetric  base  flow  and  the  other  for 
calculating  the  unsteady  three  dimensional  flow. 

4.8  Boundary  Conditions 

For  the  inflow  boundary,  the  axial  and  azimuthal  velocity  and  the  tem¬ 
perature  are  set  to  a  constant  value.  In  a  subsonic  calculation,  either  a 
Dirichlet  condition  for  pressure  and  turbulent  kinetic  energy  and  dissipation 
rate  and  a  Neumann  condition  on  radial  velocity  can  be  set  or  vice  versa. 
For  supersonic  simulations  the  pressure  is  fixed  in  the  supersonic  region  and 
a  Neumann  condition  is  used  in  the  subsonic  part  of  the  boundary  layer. 
The  outflow  boundary  is  computed  via  one-sided  stencils,  setting  the 


second  derivative  to  zero.  Only  for  subsonic  cases,  the  pressure  can  be  set  to 
a  constant  value. 
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The  freestream  boundary  is  set  according  to  Thompson  (1987)  for  super¬ 
sonic  cases.  In  subsonic  cases,  this  siplifies  to  Neumann  conditions. 

At  the  walls,  no  slip  and  no  penetration  is  enforced  on  all  velocity  com¬ 
ponents  and  the  turbulent  kinetic  energy  is  set  to  zero.  Pressure  can  be 
either  computed  at  the  wall  through  the  wall-normal  momentum  equation  or 
a  Neumann  condition  can  be  applied.  Temperature  can  be  set  to  a  constant 
value  or  computed  through  a  high-order  extrapolation,  i.e.  either  isothermal 
or  adiabatic  walls  are  specified. 

For  all  boundaries  mentioned  above,  the  density  is  computed  through  the 


equation  of  state  for  a  perfect  gas. 
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5.  CODE  VALIDATION 

5.1  Navier-Stokes  Code  for  DNS 

For  validation  of  the  first  component  of  the  Navier-Stokes  code  (for  cal¬ 
culating  the  steady  axisymmetric  base  flow),  an  axisymmetric  wake  for  a 
free  stream  Mach  number  of  M  =  0.2  and  a  global  Reynolds  number  of 
R^d  =  1, 000  was  calculated  and  compared  to  results  from  a  calculation  us¬ 
ing  an  incompressible  Navier-Stokes  code  developed  by  Schwarz  et  al.  (1994) 
[Rcd  =  {uooD)/i/oo,  where  D  is  the  diameter  of  the  body].  For  the  incom¬ 
pressible  simulations,  the  Navier-Stokes  code  of  Schwarz  et  al.  (1994)  is  based 
on  a  vector  potential  vorticity  formulation.  In  addition,  the  numerical  com¬ 
ponents  of  the  code  are  also  very  different  from  the  ones  used  in  our  present 
code.  Comparisons  between  the  compressible  and  incompressible  simula¬ 
tions  have  shown  very  good  agreement.  Comparison  of  a  global  quantity, 
the  length  of  the  recirculation  zone,  and  a  local  quantity  of  the  flow  field, 
the  azimuthal  vorticity  at  the  corner  of  the  base,  showed  practically  identical 
results.  For  validation  of  the  other  component  of  the  Navier-Stokes  code  for 
calculating  the  unsteady  (disturbed)  flow,  the  response  of  the  steady  sub¬ 
sonic  flow  field  for  M  =  0.2  and  Rcd  =  1,000  (obtained  from  the  steady 
calculation  as  discussed  above)  to  a  three  dimensional  disturbance  input  was 
calculated  and  compared  to  calculations  by  Schwarz  et  al.  (1994).  The  dis- 
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turbance  was  generated  by  a  single  pulse  at  one  location  in  the  near  wake 
region.  With  this  disturbance  generation,  a  broad  spectrum  of  disturbance 
frequencies  is  introduced  into  the  flow  field.  From  the  incompressible  simu¬ 
lations  by  Schwarz  et  al.  (1994),  it  was  found  that,  for  the  global  Reynolds 
number  of  Re  =  1, 000,  the  flow  field  includes  a  region  of  absolute  instability 
with  regard  to  three  dimensional  (helical)  disturbances.  For  additional  vali¬ 
dation  of  the  DNS  code,  we  calculated  the  response  of  the  flow  for  the  same 
conditions  as  above  (M  =  0.2,  Re^  =  1, 000)  to  strong,  nonlinear,  continuous 
excitations  through  a  blowing  and  suction  slot.  The  disturbances  introduced 
were  purely  in  the  first  helical  Fourier  mode  (k=l).  The  results  were  directly 
compared  to  experiments  that  were  conducted  in  our  hydrodynamics  labora¬ 
tory  water  channel.  A  typical  response  of  the  flow  field  is  shown  in  Figure  5.1, 
where  instantaneous  contour  lines  of  constant  azimuthal  vorticity  are  shown. 
In  this  view,  the  flow  field  exhibits  energetic  alternating  structures  (with  an 
amplitude  of  15%  of  the  free  stream  velocity)  that  are  strongly  periodic.  In 
fact.  Figure  5.2,  which  is  for  the  same  time  instant,  a  view  of  a  plane  90° 
from  that  in  Figure  5.1,  confirms  that  this  instantaneous  flow  field  consists 
of  two  counter-rotating  helical  modes.  Numerical  flow  visualization  obtained 
from  the  DNS  data,  as  shown  in  Figure  5.3  and  5.4,  was  also  compared  to 
physical  flow  visualization  of  the  laboratory  experiments.  The  agreement 
between  numerical  and  experimental  observations  was  remarkable.  As  in  the 
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Figure  5.1  Isocontours  of  constant  azimuthal  vorticity  for  M  =  0.2  and 
Re  =  1,000,  e  =  0°. 

experiments,  the  frequency  of  the  developing  structures  also  locked  in  to  that 
of  the  forcing  frequency. 

5.2  Reynolds  Averaged  Navier  Stokes 

The  turbulence  model  was  first  validated  with  a  flat  plate,  zero  pressure 
gradient  turbulent  boundary  layer  because  of  the  available  data  for  this  case 
and  the  experience  our  research  group  has  had  with  that  case  in  the  past. 
Figure  5.5  shows  the  streamwise  velocity  in  wall  coordinates  to  emphasize 
the  near-wall  behavior.  The  simulations  were  run  using  the  ASM  in  domain 
1  (see  figure  4.5)  at  M  =  0.2  with  the  incompressible  coefficients  ctj  given  in 
section  3.4.  As  initial  condition,  a  laminar  boundary  layer  profile  was  cho- 


Figure  5.2  Isocontours  of  constant  azimuthal  vorticity  for  M  —  0.2  and 
Re  =  1, 000,  e  =  90° 


Figure  5.3  Flow  visualization  using  particles.  From  DNS  data:  M  =  0.2  and 
Re  =  1,000,  e  =  0°. 
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Figure  5.4  Flow  visualization  using  particles.  From  DNS  data:  M  =  0.2  and 
Re  =  1,000,  e  =  90° 


Figure  5.5  Streamwise  velocity  in  wall-coordinates  for  turbulent  boundary 
layer  (incompressible  model). 


sen  for  the  velocity  field  and  both  turbulent  quantities  K  and  e  were  set  to 
a  constant  (typically  10“^)  throughout  the  domain.  A  derivative  condition 
was  prescribed  for  K  and  e  at  the  inflow.  This  method  was  chosen  to  ensure 
that  the  numerical  procedure  is  capable  not  only  of  maintaining  a  correct 
solution  but  also  of  producing  a  meaningful  result  from  conditions  that  differ 
significantly  from  the  desired  final  state.  The  results  of  an  ASM  case  using 
the  ramping  function  /e2  and  computations  of  the  wall  distance  independent 
ASM  (see  equation  3.48),  varying  the  constant  Ct,  are  compared  to  the  the¬ 
oretical  curves  (laminar  sublayer  with  C/+  =  and  logarithmic  law  region 
with  U'^  =  2.5log{y'^)  +  5.1)  and  a  case  run  by  D.  von  Terzi  using  a  com¬ 
pressible  code  in  Cartesian  coordinates  written  in  C  (NSCC).  All  simulations 
share  a  slope  that  is  slightly  too  steep,  but  otherwise  the  agreement  to  the 
theoretical  curves  is  satisfactory. 

5.3  ASM  Results  for  Subsonic  Wakes 

Initially,  the  turbulence  model  was  tested  in  ’’full  wake  mode”  for  an  in¬ 
compressible  case  at  M  =  0.2  with  Reo  =  100, 000  to  demonstrate  the  ability 
of  the  ASM  to  capture  unsteady  structures  (the  incompressible  case  develops 
larger  structures  due  to  the  lack  of  the  damping  due  to  compressibility).  A 
typical  result  is  shown  in  figure  5.6.  Shown  are  instantaneous  contours  of 
vorticity  and  it  is  clearly  visible  that  the  shear  layer  develops  a  strong  insta- 
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Figure  5.6  Instantaneous  contours  of  azimuthal  vorticity  for  2D-URANS; 
Ren  =  100, 000;  M  =  0.2. 


Figure  5.7  3-dimensional  distribution  of  turbulent  kinetic  energy  for  3D- 
URANS;  Reo  =  100, 000;  M  =  0.2. 

bility  with  roll-up  of  vortices  even  if  the  turbulence  model  is  fully  switched 
on  and  no  external  forcing  is  applied. 

Several  3-D  URANS  computations  were  carried  out  for  the  same  subsonic 
test  case.  Figure  5.7  clearly  shows  that  the  distribution  of  turbulent  kinetic 
energy  has  significant  variation  in  azimuthal  direction  and,  therefore,  justifies 
the  approach  of  coding  all  turbulent  variables  fully  three-dimensional. 
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5.4  Flow  Simulation  Methodology 

In  addition,  we  have  implemented  the  new  Flow  Simulation  Methodology 
(FSM)  together  with  the  new  turbulence  models  that  are  an  integral  part 
of  FSM.  We  have  tested  and  validated  almost  all  of  the  critical  individual 
modules  of  the  FSM. 

Typical  results  from  validating  the  individual  modules  of  FSM  are  pre¬ 
sented  in  Figures  5.8  and  5.9.  These  figures  were  selected  to  demonstrate 
the  role  of  the  contribution  function f  {A/ Lk),  (equation  3.50),  which  is  the 
crucial  element  of  the  FSM.  An  incompressible  turbulent  flat-plate  boundary 
layer  was  chosen  for  this  validation  because  a  wealth  of  experimental,  theo¬ 
retical,  and  numerical  data  is  available  for  comparison.  Shown  in  Figures  5.8 
and  5.9  are  color  contours  of  instantaneous  spanwise  vorticity  and  overlaid 
are  isolines  of  instantaneous  kinetic  energy.  The  results  in  Figures  5.8  and 
5.9  are  from  FSM  calculations  with  different  values  of  (equation  3.50).  A 
large  value  of  /?  leads  to  a  larger  contribution  of  the  turbulence  model  to  the 
stresses,  which  in  essence  is  equivalent  to  carrying  out  a  simulation  with  a 
coarser  grid  [for  Figure  5.9].  Therefore,  comparing  Figures  5.8  and  5.9,  it  is 
obvious  that  details  of  the  flow  due  to  smaller  structures  are  removed  when 
the  contribution  of  the  model  is  increased.  Nevertheless,  the  high  coherence 
of  vorticity  and  turbulent  kinetic  energy  is  remarkable,  even  as  fewer  details 


Streamwise  Direction,  Rejj*10’® 


Figure  5.8  FSM  applied  to  turbulent  flat-plate  boundary  layer  (incompress¬ 
ible).  Color  contours  of  instantaneous  spauiwise  vorticity  and  isolines  of  in¬ 
stantaneous  turbulent  kinetic  energy:  FSM:  ^  =  0.02. 

of  the  instantaneous  flow  field  are  resolved. 

As  discussed  in  section  3.5,  with  /(A/Lk),  the  contribution  of  the  turbu¬ 
lence  model  adjusts  locally  and  instantaneously  during  the  simulation  from 
which  Figures  5.8  and  5.9  were  obtained.  For  example,  for  the  same  bound¬ 
ary  layer  simulation.  Figure  5.10  shows  the  time-averaged  /{A/Lk)-  It  is 
obvious  that  the  contribution  function  behaves  as  expected.  The  contribu¬ 
tion  of  the  model  is  large  close  to  the  wall,  where  the  physical  resolution 
(grid  size  compared  to  Kolmogorov  length  scale)  is  coarse.  The  contribution 
then  decreases  since  the  Kolmogorov  scale  increases  with  increasing  distance 
from  the  wall.  This  is  a  clear  demonstration  that  the  contribution  function, 
equation  5.10,  is  performing  as  it  was  designed  to  perform. 


Wall-Normal,  y/L 


Streamwise  Direction,  Rejj*10’® 


Figure  5.9  FSM  applied  to  turbulent  flat-plate  boundary  layer  (incompress¬ 
ible).  Color  contours  of  instantaneous  spanwise  vorticity  and  isolines  of  in¬ 
stantaneous  turbulent  kinetic  energy:  FSM:  ^  =  0.05. 


Figure  5.10  Time-averaged  profiles  of  /(A/L^)  relative  to  local  boundary 
layer  thickness  for  cases  shown  in  Figures  5.8  and  5.9  and  for  larger  values 
of/(A/LK). 
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6.  RESULTS 

6.1  DNS 

6.1.1  Steady  (Undisturbed  Flow) 

We  have  performed  supersonic  calculations  for  a  Mach  number  of  M  = 
2.46  [experiments  of  Herrin  &  Dutton  (1991)  and  the  Reynolds  averaged 
calculations  by  Sahu  (1992)].  We  made  calculations  for  different  Reynolds 
numbers  in  order  to  explore  the  effect  of  Reynolds  number  on  the  resulting 
flow  field.  For  the  same  Reynolds  number  as  for  the  subsonic  case,  Reo  = 
1,000,  the  base  pressure  was  found  to  be  much  lower  for  the  supersonic 
case  than  for  the  subsonic  case,  which  is  consistent  with  the  experimental 
observations  of  Herrin  k,  Dutton  (1991).  However,  in  the  experiments,  the 
pressure  distribution  was  found  to  be  practically  constant  along  the  base, 
while  it  varied  by  almost  20%  in  the  numerical  results.  A  variation  of  pressure 
along  the  base  was  also  found  in  the  Reynolds  averaged  calculations  (using 
turbulence  models)  by  Sahu  (1992)  and  by  Sturek  k  Nietubicz  (1992).  One 
may  speculate  that  the  different  behavior  of  the  base  pressure  along  the  base 
is  due  to  the  action  of  the  large  structures,  which  is  not  included  in  our 
steady  base  flow  calculations  and  may  not  be  adequately  represented  in  the 
turbulence  models  used  in  the  Reynolds  averaged  calculations.  In  fact,  results 
from  our  numerical  simulations  discussed  below  support  this  conjecture.  An 
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Pressure 


Figure  6.1  Contours  of  constant  pressure  for  Re  =  18,000,  M  =  2.46. 

answer  to  this  speculation  is  of  considerable  relevance  for  the  development 
of  improved  turbulence  models.  In  order  to  obtain  a  closer  comparison  with 
the  experimental  data  by  Herrin  &:  Dutton  (1991),  we  also  calculated  steady 
flow  fields  for  Reynolds  numbers  up  to  Ren  =  100,000,  with  a  boundary 
layer  thickness  at  the  corner  of  5  =  0.07D,  which  matched  the  boundary 
layer  thickness  of  the  experiment.  A  typical  result  for  Ren  =  18, 000  in  the 
form  of  constant  pressure  contours  is  shown  in  Figure  6.1. 

6.1.2  Unsteady  (Disturbed)  Three-Dimensional  Flow. 

To  investigate  the  absolute  stability  behavior  of  the  supersonic  base  flow, 
a  pulse  disturbance  was  introduced  into  the  recirculation  region.  For  the  first 
calculation,  we  used  the  same  Reynolds  number  as  for  the  incompressible  cal¬ 
culation  {  Rsd  =  1,000).  After  an  initially  very  short  strong  response  of  the 
flow  field,  the  disturbances  decayed  everywhere.  Increasing  the  Reynolds 
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number,  2, 000  <  Reo  <  25, 000,  showed  no  difference  in  the  disturbance  re¬ 
sponse:  disturbances  decayed  after  a  relatively  short  time  with  the  exception 
that  the  decay  rate  decreased  with  increasing  Reynolds  number.  Thus,  in 
contrast  to  the  subsonic  case,  where  an  absolute  instability  exists  with  respect 
to  helical  disturbance  modes,  the  calculations  for  Mqo  =  2.46  and  Reynolds 
numbers  up  to  Reo  =  25,000  showed  no  indication  of  absolute  instability 
with  respect  to  helical  or  axisymmetric  modes.  Therefore,  it  appears,  as 
has  been  expected,  that  supersonic  wake  flows  are  more  stable  than  their 
incompressible  counterparts.  This  conjecture  is  consistent  with  the  findings 
of  Chen  et  al.  (1990)  for  a  two  dimensional  compressible  wake,  where  the 
growth  rates  of  disturbances  were  found  to  be  up  to  an  order  of  magnitude 
smaller  than  for  a  comparable  incompressible  wake  and  structures  appeared 
to  diffuse  much  faster. 

6.1.3  Search  for  Absolute  Instability  for  Supersonic  Axisym¬ 
metric  Base  Flows 

As  seen  from  the  previous  discussion,  the  results  for  the  supersonic  flow 
field  of  M  =  2.46  and  Reynolds  numbers  up  to  Re^  =  25, 000  do  not  exhibit 
an  absolute  instability.  Since  compressibility  has  a  stabilizing  effect  (as  dis¬ 
cussed  above),  we  decided  to  perform  a  simulation  for  a  lower  Mach  number 
to  see  if  an  absolute  instability  would  arise.  Thus,  we  calculated  the  flow 
field  for  a  Mach  number  of  M  =  1.2  and  Re^  =  4,000  (Tourbier  &  Fasel 
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Figure  6.2  Time  response  of  radial  momentum  at  z  =  0.5  and  r  =  0.5 
for  three  different  azimuthal  Fourier  modes  (k)  to  a  three-dimensional  pulse 
disturbance  in  the  near  wake  for  M  =  1.2  and  Rcd  =4, 000. 

(1994)).  As  before,  the  flow  field  was  disturbed  by  introducing  a  single  pulse 
locally  in  the  recirculation  region.  A  typical  temporal  response  of  the  flow 
field  is  given  in  Figure  6.2.  It  is  obvious  that,  for  different  azimuthal  Fourier 
modes,  the  disturbance  grows  exponentially  in  time  until  it  reaches  a  state 
of  nonlinear  saturation.  This  is  a  clear  manifestation  of  the  existence  of  an 
absolute  instability. 

An  impression  of  the  instantaneous  flow  field  can  be  obtained  from  Figures 
6.3  and  6.4.  Shown  are  instantaneous  contour  lines  of  the  azimuthal  vorticity 
after  the  disturbance  has  reached  the  state  of  nonlinear  saturation.  Large 


structures  can  be  observed  in  the  near-wake  region.  In  order  to  investigate 
the  influence  of  the  structures  on  the  global  flow  field,  the  (time-averaged) 
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Figure  6.3  Isocontours  of  constant  azimuthal  vorticity  for  M  —  1.2  and 
i?e  =  4,000,  0  =  0°. 


-1.0  -0.5  o!o  0^5  l!o  1.5  2!o  2!5  3.0  3.5  4^0  4.5 

X 


Figure  6.4  Isocontours  of  constant  azimuthal  vorticity  for  M  =  1.2  and 
Re  =  4, 000,  e  =  90° 
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Figure  6.5  Base  pressure  distribution  for  a  steady  axisymmetric  flow  field 
and  a  time-averaged  unsteady  flow  field  for  M  =  1.2  and  Re  =  4, 000. 

mean  flow  was  calculated.  A  comparison  of  the  base  pressure  obtained  from 
a  steady  calculation  and  that  obtained  from  a  time-averaged  mean  flow  cal¬ 
culation  is  shown  in  Figure  6.5.  Clearly,  the  presence  of  the  large  structures 
causes  a  drop  in  the  base  pressure.  Also,  for  the  unsteady  calculation,  the 
distribution  of  the  pressure  along  the  base  is  virtually  constant,  which  agrees 
with  observations  in  experiments,  while  pressure  varies  strongly  for  the  steady 
calculation.  These  results  are  a  strong  indication  that  the  action  of  the  large 
structures  is  responsible  for  the  flat  base  pressure  distribution  observed  in 
experiments.  Thus,  there  is  further  evidence  that  this  may  also  be  the  reason 
why  Reynolds-averaged  calculations  have  difficulties  in  reproducing  the  flat 
pressure  distribution. 

In  other  simulations  for  M  =  2.46,  we  pushed  the  Reynolds  number 
to  Rer)  =  100,000.  As  before,  the  flow  field  was  disturbed  by  a  single 
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Figure  6.6  Isocontours  of  constant  azimuthal  vorticity  for  M  =  2.46  and 
i2e  =  30,000,  0  =  0°. 

pulse  in  the  recirculation  region.  Now,  as  for  the  M  =  1.2  case  and  in 
contrast  to  the  case  with  M  =  2.46  and  Rcd  =  25,000,  the  disturbances 
grow  exponentially  in  time  (starting  at  i?e£)  =  30, 000),  a  clear  indication  of 
an  absolute  instability.  Figures  6.6  and  6.7  show  instantaneous  contour  lines 
of  total  vorticity  for  Re^  —  30, 000.  It  is  obvious  that  the  large  structures 
exist;  in  fact,  we  found  that  the  intensity  (amplitude)  of  these  structures  is 
on  the  order  of  15%  of  the  free  stream  velocity.  Thus,  they  are  at  least  of 
the  same  relative  strength  as  for  incompressible  wakes. 

The  qualitative  effect  of  the  dynamics  of  the  large  structures  on  the  time- 
averaged  base  pressure  is  evident  from  Figure  6.8.  As  for  the  lower  Mach 
number  case  (M  =  1.2;  see  above),  the  base  pressure  of  the  time-averaged 
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Figure  6.7  Isocontours  of  constant  azimuthal  vorticity  for  M  —  2.46  and 
Re  =  30, 000,  d  =  90° 

flow  field  has  an  almost  constant  distribution  over  the  radius,  which  is  con¬ 
sistent  with  experimental  results.  In  this  case,  the  effect  of  base  pressure 
reduction  due  to  the  time-dependent  structures  is  even  larger  than  for  the 
lower  Mach  number  case. 

For  Reynolds  number  larger  than  30,000  the  old  fourth-order  explicit 
code  required  exceedingly  costly  computations,  so  all  DNS  conducted  for 
Reu  >  30, 000  employed  the  new,  higher  order  code.  Figure  6.9  shows  color- 
contours  of  total  vorticity  of  the  plane  ^  =  0°  for  a  full  3-D  calculation  at 
Rei)  =  100,000.  One  can  clearly  see  that  the  simulation  employing  the  new 
sixth-order  compact  difference  stencils  is  able  to  resolve  small  scale  structures 
throughout  the  domain  shown.  Also,  it  is  important  to  note  that  the  flow 
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Figure  6.8  Base  pressure  distribution  for  a  steady  axisymmetric  flow  field 
and  a  time-averaged  unsteady  flow  field  for  M  =  1.2  and  Re  =  4, 000. 

behavior  at  the  points  on  the  axis-line  do  not  appear  explicitly  in  this  picture 
but  rather  are  like  at  any  other  point  in  the  field,  an  indication  that  the  new 
axis-treatment  is  working  very  well. 

End-view  images  of  instantaneous  total  vorticity  are  displayed  in  figure 
6.10  for  different  downstream  locations.  Figure  6.10  a)  is  a  slice  very  close 
to  the  base  and  shows  the  small  structures  in  the  recirculation  zone  and  no 
apparent  symmetry  is  present  (except  that  with  respect  to  0  =  0°).  At  the 
locations  b)  and  c)  which  are  comparable  to  locations  C  and  D  in  Bour¬ 
don  &  Dutton  (1998),  turbulent  structures  that  have  undergone  substantial 
growth  can  be  seen.  Also,  the  presence  of  longitudinal  structures  in  the  shear 
layer  is  clearly  visible.  The  position  d),  which  is  a  cut  in  the  trailing  wake 
and  can  be  compared  to  position  E  in  Bourdon  &  Dutton  (1998),  shows  the 
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Figure  6.9  Instantaneous  color-contours  of  total  vorticity  for  3D-DNS;  Reo  = 
100,000;  M  =  2.46;  0  =  0°. 

formation  of  a  symmetry  with  respect  to  0  =  90°  as  observed  in  the  experi¬ 
ments.  This  clearly  indicates  that  even  though  the  Reynolds  number  chosen 
for  the  simulations  is  an  order  of  magnitude  smaller  than  that  in  the  exper¬ 
iment,  qualitative  agreement  is  achieved.  Note  that  all  results  shown  here 
are  instantaneous  views  of  the  unsteady  flow  field. 

Figures  6.11  and  6.12  show  instantaneous  contour  lines  of  total  vorticity 
for  Reo  =  100,000.  Again,  it  is  visible  that  large  structures  exist  and  it  is 
clear  that  the  higher  order  code  resolves  more  scales  than  the  formerly  used 
fourth-order  explicit  code.  This  is  also  impressively  demonstrated  in  figures 
6.13  and  6.14  where  instantaneous  contours  and  colour-contours  are  shown 
for  Mach  lines  and  pressure,  respectively.  The  new  code  is  able  to  capture 
the  shocklets  that  occur  in  the  flow,  originating  from  the  structures  in  the 
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c)  z=8  d)  z=12 

Figure  6.10  Contours  of  instantaneous  total  vorticity  for  various  downstream 
locations;  Rev  =  100,000;  M  =  2.46. 
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Figure  6.11  Isocontours  of  constant  azimuthal  vorticity  for  M  =  2.46  and 
i2e  =  100,000,  0  =  0°. 


Figure  6.12  Isocontours  of  constant  azimuthal  vorticity  for  M  =  2.46  and 
iile  =  100, 000,  0  =  90°. 

shear  layer  and  travelling  downstream.  Resolving  these  pressure-gradients  is 
necessary  for  describing  the  dynamics  of  the  flow  correctly. 

Another  proof  of  the  existence  of  large  structures  was  also  established  in 
our  Direct  Numerical  Simulations  of  a  two-dimensional  base  flow  for  M  — 
2.46  and  Rbd  =  250, 000.  Typical  results  from  these  investigations  are  shown 
in  Figures  6.15  and  6.16.  It  is  obvious  that  large  structures  are  present,  even 
when,  as  in  this  case  (larger  Reynolds  number),  they  are  resulting  from  an 
absolute  instability  and  did  not  require  continuous  forcing.  The  structures 
are  most  dynamic  in  the  shear  layer  (see  Figures  6.15  and  6.16),  where  vortex 
merging  occurs  and  results  in  larger,  more  energetic  structures. 


90 


density 


Figure  6.15  Contours  of  instantaneous  density  for  DNS  of  two-dimensional 
base  flow  for  M  =  2.46  and  Re  =  250, 000. 


Figure  6.16  Contours  of  instantaneous  spanwise  vorticity  for  DNS  of  two- 
dimensional  base  flow  for  M  =  2.46  and  Re  =  250, 000. 
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6.2  Large  Eddy  Simulation 

For  investigations  of  higher  Reynolds  number  turbulent  flows  and  as  a 
benchmark  for  unsteady  RANS  and  FSM  calculations,  we  also  implemented 
subgrid-scale  turbulence  models  for  Large-Eddy  Simulations  (LES).  In  a  first 
step,  we  implemented  a  Smagorinsky-type  subgrid-scale  turbulence  model 
into  our  original  DNS  code.  This  turbulence  model  was  suggested  for  com¬ 
pressible  flow  by  Professor  C.  Speziale  of  Boston  University.  With  this  model, 
we  calculated  turbulent  wake  flows  for  a  range  of  Reynolds  numbers,  from 
transitional  all  the  way  to  fully  turbulent  wakes. 

In  order  to  test  and  validate  the  LES  code,  we  have  performed  simulations 
for  the  same  Mach  number,  M  =  2.46,  as  in  the  experiments  carried  out  at 
the  University  of  Illinois  (Dutton  and  co-workers).  However,  the  range  of 
Reynolds  numbers  in  our  simulations  was  considerably  lower,  between  30,000 
and  400,000  (based  on  base  diameter).  Typical  results  of  LES  obtained 
with  the  standard  fourth-order-accurate  code  and  with  the  Smagorinsky- 
type  subgrid-scale  model  (as  proposed  for  compressible  flow  by  C.  Speziale) 
are  presented  in  Figures  6.17  -  6.32  for  M  =  2.46  and  Reynolds  numbers 
Rcd  =  30, 000  and  100,000. 

Comparison  of  corresponding  plots  for  Reynolds  numbers  30,000  and 
100,000  (Figures  6.17  -  6.32)  indicates  that  there  are  still  quantitative  and 
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Figure  6.17  Isolines  of  instantaneous  total  vorticity  from  LES  with  a 
Smagorinsky  model:  M  =  2.46  and  Re  =  30, 000. 
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Figure  6.18  Isolines  of  instantaneous  total  vorticity  from  LES  with  a 
Smagorinsky  model:  M  =  2.46  and  Re  =  100, 000. 

qualitative  changes  of  the  flow  fleld  and  the  turbulent  statistics  as  the  Reynolds 
number  increases  from  30,000  to  100,000.  However,  when  comparing  the  re¬ 
sults  for  Reynolds  numbers  between  100,000  and  400,000  (not  shown  here), 
the  changes  with  Reynolds  number  are  not  significant,  which  may  be  an  indi¬ 
cation  that  the  accuracy  of  the  code  together  with  the  Smagorinsky  subgrid- 
scale  model  was  not  sufficient  for  high-quality  simulations  that  would  allow 
detailed  investigation  of  the  flow  physics  for  high  Reynolds  numbers.  This 
was  the  motivation  for  developing  the  considerably  more  accurate  compact 
codes  and  for  developing  the  Flow  Simulation  Methodology  (FSM)  (as  dis¬ 
cussed  in  3.5)  so  that  reliable  simulations  can  be  performed  for  considerably 
higher  Reynolds  numbers. 
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Figure  6.19  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  30, 000  from  LES  with  a  Smagorinsky  model;  pressure. 


Figure  6.20  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  30, 000  from  LES  with  a  Smagorinsky  model:  density. 


Figure  6.21  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  30, 000  from  LES  with  a  Smagorinsky  model:  radial  velocity. 
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Figure  6.22  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
30, 000  from  LES:  radial  turbulence  intensity. 


Figure  6.23  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
30, 000  from  LES:  azimuthal  turbulence  intensity. 


Figure  6.24  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
30, 000  from  LES:  azimuthal  component  of  turbulent  Reynolds  stress. 


Figure  6.25  Isolines,  for  time-averaged  quantities  with  M  =  2.46  and  Re 
30, 000  from  LES:  turbulent  kinetic  energy. 
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Figure  6.26  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  100, 000  from  LES  with  a  Smagorinsky  model:  pressure. 


Figure  6.27  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  100,  000  from  LES  with  a  Smagorinsky  model:  density. 


Figure  6.28  Isolines  for  the  time-averaged  flow  field  with  M  =  2.46  and 
Re  =  100, 000  from  LES  with  a  Smagorinsky  model:  radial  velocity. 


Figure  6.29  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
100, 000  from  LES:  radial  turbulence  intensity. 


Figure  6.30  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
100, 000  from  LES:  azimuthal  turbulence  intensity. 


k> 


°ao 


-0.001 


T!o  Tls  ?o  2!$  alo  ilo  ils  So  Ss  So  Si  rlo  yis 

X 


Figure  6.31  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
100, 000  from  LES:  azimuthal  component  of  turbulent  Reynolds  stress. 
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Figure  6.32  Isolines  for  time-averaged  quantities  with  M  =  2.46  and  Re 
100, 000  from  LES:  turbulent  kinetic  energy. 
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6.3  Reynolds  Averaged  Navier  Stokes 

So  far,  only  2-D  simulations  have  been  performed  for  the  supersonic  case 
by  Dutton  and  co-workers.  The  initial  condition  for  the  supersonic,  high 
Reynolds  number  cases  was  the  axisymmetric  result  for  Reo  =  100, 000  and 
M  =  2.46  obtained  by  DNS.  In  addition,  the  distribution  of  K  and  e  was  set 
to  a  uniform  constant  throughout  the  computational  domain.  All  calculations 
were  performed  on  a  relatively  coarse  grid  with  250  points  in  the  streamwise 
and  120  points  in  the  radial  direction. 

Figures  6.33-6.36  show  the  influence  of  various  inflow  boundary  conditions 
for  K.  When  K  is  kept  constant  at  the  inflow  (6.33),  the  turning  angle 
of  the  shear  layer  is  larger  and  the  recirculation  length  shorter  than  if  a 
Neumann  condition  is  used  (cases  6.34  and  6.35).  This  can  be  explained  by 
the  fact  that  for  the  ASM  model,  K  and  e  vanish  in  the  approach  boundary 
layer  if  K  is  not  kept  at  a  constant  value  at  the  inflow,  thus  producing  a 
laminar  boundary  layer  that  generates  unsteady  shedding  in  the  shear  layer. 
In  case  6.36,  the  standard  K  —  s  model  produces  a  more  turbulent  boundary 
layer  and,  consequently,  comes  closest  to  the  experiments  by  Dutton  and 
co-workers  in  terms  of  turning  angle  and  recirculation  length. 


-1.0  -2.0  2.0  1.0  0.0  -1.0  -2.0 


98 


mrmm 


-4.0  -2.0  0.0  2.0  4.0 


Vorticity  Magnitude 


-1.0  0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0  10.0 

z 


Figure  6.33  Color-contours  of  instantaneous  vorticity  for  2-D  ASM,  with 
Dirichlet  condition  for  K  at  inflow;  Rcd  =  3, 300, 000;  M  =  2.46. 
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Figure  6.34  Color-contours  of  instantaneous  vorticity  for  2-D  ASM,  with 
Neumann  condition  for  K  at  inflow;  Rcd  =  3,300,000;  M  =  2.46. 
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Figure  6.35  Color-contours  of  averaged  vorticity  for  2-D  ASM,  with  Neumann 
condition  for  K  at  inflow  ;  Reo  =  3,300,000;  M  =  2.46. 
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Figure  6.36  Color-contours  of  instantaneous  vorticity  for  2-D  standard  K-e- 
model;  Rcd  =  3,300,000;  M  =  2.46. 
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6.4  Flow  Simulation  Methodology 

Finally,  FSM  has  been  applied  to  the  wake-flow,  even  though  so  far  only 
in  2-D.  Figures  6.37  -  6.40  show  color  contours  of  instantaneous  vorticity 
for  a  test-case  at  Reo  —  100,000  and  M  =  0.2  for  a  2-D  ”DNS”,  2-D 
”FSM”  with  different  P  and  a  2-D  ’’URANS”  calculation,  all  with  a  resolution 
210x120x1.  It  is  obvious  that  the  amount  of  dissipation  becomes  larger  for 
higher  contributions  of  the  Reynolds  Stresses.  In  fig.6.37,  no  turbulence 
model  is  used  and  because  of  the  simulation  being  2-dimensional,  there  is 
a  significant  amount  of  energy  trapped  in  the  2-D  mode  downstream  of  the 
corner. 

The  strength  of  that  vortex  even  deforms  the  shear  layer  and  generates 
a  massive  shear  layer  break-up  close  to  the  base  {z  1.5).  For  the  FSM 
calculation  with  /?  =  0.0001  shown  in  fig.  6.38,  the  instability  in  the  shear 
layer  is  damped  due  to  the  higher  amount  of  dissipation  and  the  vorticies 
develop  further  downstream  (z  «  2.5). 

The  contribution  function  in  this  case  gave  values  between  5-25%.  For  the 
case  of  /?  =  0.002  and  the  URANS  calculation  shown  in  figures  6.39  and  6.40, 
respectively,  the  instantaneous  flow  field  looks  very  similar,  largely  because 
for  this  value  of  (3  and  any  larger  value  the  contribution  of  the  turbulence 
model  is  very  large  in  the  whole  domain  and  therefore  the  RANS-limit  is 
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Figure  6.37  Color-contours  of  instantaneous  vorticity  for  2-D  ”DNS”;  Reo  = 
100,000;  M  =  0.2. 

approached.  In  these  cases  the  vorticies  develop  even  further  downstream 
(z  «  4.0)  and  there  is  no  deformation  of  the  shear  layer  anymore.  However, 
it  is  important  to  note  that  the  simulation  still  produces  these  unsteady 
structures  even  with  the  large  amount  of  turbulent  dissipation  added  (no 
forcing  is  required  to  obtain  these  structures;  they  arise  due  to  the  strong 
shear  layer  instability). 


Figure  6.38  Color-contours  of  instantaneous  vorticity  for  2-D  ”FSM”  with 
P  =  0.0001;  Rbd  =  100,000;  M  =  0.2. 


Figure  6.40  Color-contours  of  instantaneous  vorticity  for  2-D  URANS; 
ReD  =  100,000;  M  =  0.2. 
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7.  CONCLUSIONS 

For  bluff  bodies  at  supersonic  speeds,  a  significant  portion  of  the  drag 
is  associated  with  base  drag,  caused  by  the  low  pressure  on  the  base.  This 
base  pressure  is  greatly  influenced  by  large  vortical  structures,  which,  for 
axisymmetric  bodies,  are  of  helical  nature.  While  the  vortex  formation  pro¬ 
cess  and  the  nature  of  the  helical  structures  has  been  studied  in  depth  for 
incompressible  flows,  both  experimentally  and  numerically,  considerably  less 
is  known  regarding  their  information  for  supersonic  flows.  Experiments  have 
confirmed  the  existence  of  large  structures  in  supersonic  flows,  but  their  be¬ 
havior  is  poorly  understood. 

A  new,  high-order  compressible  Navier-Stokes  code  using  compact  dif¬ 
ference  stencils  of  6*^-order  accuracy  derived  for  non-equidistant  grids  and 
a  state-of-the-art  axis  treatment  was  developed  to  investigate  the  stability 
behavior  in  an  axisymmetric,  bluff  body  wake  at  M  =  2.46.  This  Mach 
number  was  chosen  to  enable  comparison  with  experiments  conducted  at  the 
University  of  Illinois  at  Urbana  Champaign  by  Dutton  and  co-workers. 

Turbulence  models,  such  as  Large  Eddy  Simulation  (LES)  and  Reynolds 
Averaged  Navier  Stokes  (RANS)  calculations  were  also  incorporated  into  the 
code.  The  K  —  e  equations  were  implemented  in  fully  three-dimensional 
form  to  ensure  that  RANS  simulations  can  capture  azimuthal  variations  in 
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the  flow  fleld  (e.g.,  helical  structures). 

In  particular  we  have  also  developed  and  successfully  implemented  the 
new  Flow  Simulation  Methodology  (FSM).  The  centerpiece  of  FSM  is  a  strat¬ 
egy  to  provide  the  proper  amount  of  modeling  of  the  subgrid  scales.  This  is 
accomplished  by  a  ”  contribution  function”  which  locally  and  instantaneously 
compares  the  smallest  relevant  scales  to  the  local  grid  size.  The  contribution 
function  is  designed  such  that  it  provides  no  modeling  if  the  computation  is 
locally  well  resolved  so  that  the  computation  approaches  a  Direct  Numerical 
Simulation  (DNS)  in  the  flne  grid  limit,  or  provides  modeling  of  all  scales  in 
the  coarse  grid  limit  and  thus  approaches  an  unsteady  RANS  (URANS)  cal¬ 
culation.  In  between  these  resolution  limits,  the  contribution  function  adjusts 
the  necessary  modeling  for  the  unresolved  scales  while  the  larger  (resolved) 
scales  are  computed  as  in  traditional  Large  Eddy  Simulations  (LES). 

Our  calculations  have  shown  that  an  absolute  instability  at  M  =  1.2 
exists  already  for  Reo  =  4,000,  however,  for  M  =  2.46,  the  onset  of  the 
absolute  instability  was  not  observed  until  Rcd  >  30,000.  Preliminary  un¬ 
steady  RANS  simulations  (URANS)  for  Reo  =  3, 300, 000  at  M  =  2.46  were 
performed  and  showed  coherent  structures. 

All  the  results  obtained  so  far  have  shown  that  the  developed  Navier 
Stokes  code  functions  properly  and  that  the  high-order  compact  differencing 
is  of  great  advantage  for  supersonic  base  flow  simulations.  Furthermore,  cal- 


105 


culations  with  the  implemented  Flow  Simulation  Methodology  (FSM)  have 
demonstrated  that  this  approach  has  great  promise  for  allowing  simulations 
of  supersonic  base  flows  for  much  larger  Reynolds  numbers  than  possible  with 
conventional  Large  Eddy  Simulations. 
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